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Preface 


This document describes version 2.2 of the Field Encapsulation Library (FEL), a li- 
brary of mesh and field classes. FEL is a library for programmers — it is a “building 
block” enabling the rapid development of applications by a user. Since FEL is a library 
intended for code development, it is essential that enough technical detail be provided 
so that one can make full use of the code. Providing such detail requires some assump- 
tions with respect to the reader’s familiarity with the library implementation language, 
C++, particularly C++ with templates. We have done our best to make the explanations 
accessible to those who may not be completely C++ literate. Nevertheless, familiarity 
with the language will certainly help one’s understanding of how and why things work 
the way they do. One consolation is that the level of understanding essential for using 
the library is significantly less than the level that one should have in order to modify or 
extend the library. 

One more remark on C++ templates: Templates have been a source of both joy and 
frustration for us. The frustration stems from the lack of mature or complete implemen- 
tations that one has to work with. Template problems rear their ugly head particularly 
when porting. When porting C code, successfully compiling to a set of object files typ- 
ically means that one is almost done. With templated C++ and the current state of the 
compilers and linkers, generating the object files is often only the beginning of the fun. 
On the other hand, templates are quite powerful. Used judiciously, templates enable 
more succinct designs and more efficient code. Templates also help with code mainte- 
nance. Designers can avoid creating objects that are the same in many respects, but not 
exactly the same. For example, FEL fields are templated by node type, thus the code 
for scalar fields and vector fields is shared. Furthermore, node type templating allows 
the library user to instantiate fields with data types not provided by the FEL authors. 
This type of flexibility would be difficult to offer without the support of the language. 

For users who may be having template-related problems, we offer the consola- 
tion that support for C++ templates is destined to improve with time. Efforts such as 
the Standard Template Library (STL) will inevitably drive vendors to provide more 
thorough, optimized tools for template code development. Furthermore, the benefits 
will become harder to resist for those who currently subscribe to the least-common- 
denominator “code it all in C” strategy. 

May FEL bring you both increased productivity and aesthetic satisfaction. 
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Chapter 1 

Introduction 


The Field Encapsulation Library (FEL) is a library for representing fields and meshes. 
The fields may represent the results from simulations or experimental observations. 
The main goals of FEL are to provide: 

• a horizontal product supporting the quick development of new applications 

• an extensible framework supporting a variety of mesh and field types 

• a high-performance mesh-independent interface 

• a library supporting work with large data sets 

“Horizontal products” are reusable libraries that make the development of “vertical 
products” — applications — easier. FEL is a set of C++ classes designed to support 
not just the field and mesh types that are currently implemented, but also new mesh 
and field types in the future. A key part of the design is the development of a mesh- 
independent interface, i.e., an interface that allows users to write a single algorithm that 
will work with a variety of mesh and field types. The goal is to be able to introduce new 
types in the future and reuse algorithms written in terms of the FEL interface without 
modification. Finally, there are a number of features in FEL that are designed to make 
working with large data sets feasible for more users. 

To develop an application using FEL, you should not find it necessary to read this 
entire document. The FEL release includes a primer which provides an overview of 
many features of FEL along with sample programs. The primer should appeal to those 
who prefer to “dive in” and learn by writing actual applications. On the other hand, the 
primer, in the name of brevity, focuses on library essentials. To get a more complete 
view of the features and design philosophy of FEL, you are encouraged to read on. 

To convey the flavor of programming with FEL, we will now walk through a small 
example. 
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1.1 An FEL example 

The following program reads in a mesh and field, then finds the cell with the greatest 
pressure gradient magnitude at its centroid: 

# include "FEL.h" 

int main (int argc, char* argv[]} 

{ 

unsigned flags = FEL_deduce_mesh_type (argv[l] ) ; 

FEL_mesh_ptr mesh = FEL_read_mesh(argv[l] , flags); 
FEL_plot3d_q_f ield_ptr q_field = 

FEL_read_q {mesh, argv[2] , flags); 

FEL_f loat_f ield_ptr p_field = 

FEL_plot3d_make _pressure_f ield (q_f ield) ; 

FEL_vector3f_f ield_ptr grad_p_field = 

new FEL_gradient_of_f loat_f ieldl {p — field) ; 

FEL_f loat_f ield__ptr mag_grad_p_f ield = 

new FEL_magnitude_of_vector3f_f ield (grad_p_f ield) ; 

n\ag_grad p_f ield->set (FEL__SIMPLICIAL_DECOMPOSITION, 0) ; 

mag grad p field->set { FEL_ INTERPOLATION, 

FEL_ISOPARAMETRIC_INTERPOLATION) ; 

int res; 

float field_value, max_value = 0.0; 

FEL_vector3f coords[8]; 

FEL_phys_pos centroid; 

FEL_cell max_cell ; 

FEL_cell_i ter iter; 

for (mag grad p f ield->begin (kiter) ; ! iter . done () ;. + + iter) { 
mesh->coordinates_at_cell ( *iter , coords) ; 
centroid. set {0 . 0, 0.0, 0.0); 

for (int i = 0; i < { *iter ) . get_n_nodes ( ) ; i++) 
centroid += coords [i] ; 
centroid /= { * i ter) . get_jn_nodes { ) ; 

res = mag_grad_p_field->at_phys_pos (centroid, &f ield_value) ; 
if (res FEL_OK) continue; 

cout << "Pressure gradient magnitude at " << *iter 
<< " is " << field_value « endl; 
if (field_value > max_value) { 
max_value = field — value; 
max_cell = *iter; 

} 

cout « "Maximum of " << max_value « " in " << max_cell « endl; 
return 0; 



1.1. AN FEL EXAMPLE 


13 


The example demonstrates several key features of FEL. The first is mesh inden- 
pendence — the same program works for fields based on a structured, unstructured, or 
multi-zone mesh. The mesh independence is achieved via FEL iterators and file read- 
ing functions that construct meshes and fields of the appropriate type. The particular 
iterator used in the example, a ceil iterator, loops over the (hexahedral or tetrahedral) 
cells in a mesh. Using each cell, one can easily access mesh and field values. The 
example also illustrates how FEL makes it easy for the user to compose new fields in 
terms of existing ones. The field data in the example is known “solution data, in other 
words, data produced by a flow solver. Starting with a solution data field (cj-f iold), 
one can construct a pressure field, a gradient field and a gradient magnitude field, each 
in one statement. The program also shows how one can access field values at arbitrary 
positions using just one function call (at_phys_pos). 

There are many other capabilities beyond those shown in the example, such as 
support for time-varying fields and transformed fields. Rather than inventory ail the 
capabilities of FEL here, the remainder of this section walks through the first example 
and mentions the chapters in this document in which to look for more information. 

unsigned flags = FEL_deduce__mesh_type (argv [ 1] ) ; 

FEL_mesh__ptr mesh = FEL_read_mesh (argv [1] , flags); 

FEL_plot3d_q_f ield_ptr q_field = 

FEL_read_q (mesh, argv[2], flags); 

The first four lines read in mesh and field data from the files named by argv [ 1 ] 
and argv [ 2 ] . The program first uses the mesh type deducer function to deduce 
flags indicating the particular file type (Chapter 22). FEL currently supports reading 
PLOT3D files and “Enterprise” paged files. Mesh and field objects based on paged 
files do not load the whole data set into memory, but instead work in a demand-driven 
style where only subsets of the data are loaded as needed (Chapter 23). The user can 
also read in data from a file in some other format and then construct meshes or fields 
directly. (Chapters 12, 18). 

Fields and meshes in FEL are reference counted. The FEL type names with the 
_pt r suffix are all “smart pointers” that refer to reference counted objects (Chapter 7). 
For the most part, the pointers behave just like C-style pointers. FEL also uses C++ 
templates, which support parameterized types. For instance, fields in FEL are parame- 
terized by their node type (Chapter 17). FEL uses typedefs to hide many of the common 
template instantiations, so that most FEL programs can be written without ever using 
template syntax directly (Chapter 2). 

FEL_float_field_ptr p_field = 

FEL_plot3d_make_pressure_f ield (q_f ield) ; 

FEL_vector3f_f ield_ptr grad_p__f ield = 

new FEL_gradient__of_f loat_f ieldl (p__f ield) ; 

FEL_f loat_field_ptr mag_grad_p_f ield = 

new FEL_magnitude_of_vector3f_f ield (grad_p_f ield) ; 
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Here we create three derived fields from the fundamental solution data. The first 
field uses a built-in function to convert PLOT3D solution vectors to pressure values 
(Chapter 24). The second field applies a differential operator (“grad”) to the pressure 
field, producing a vector field from a scalar field (Chapter 20). The third field takes 
the norm of the gradient field, yielding a scalar field (Chapter 19). This progression 
illustrates how FEL allows new fields to be defined in terms of preexisting ones. 

The derived and differential operator fields require almost no storage, since all 
values are computed on demand. This demand-driven or “lazy” evaluation approach 
provides a significant advantage when working with large data sets, since the memory 
and computational requirements of precomputing a derived value over a whole field 
can be prohibitive. 

mag_grad_p_f ield->set (FEL__SIMPLICIAL_DECOMPOSITION, 0) ; 

mag_grad_p_f ield->set (FEL__INTERPOLATION, 

FEL_ISOPARAMETRIC_INTERPOLATION) ; 

These directives tell FEL how to conduct some of its numerical business (Chap- 
ters 9, 11). FEL supports several spatial interpolation modes (Chapter 9). The sim- 
plicial decomposition flag controls whether FEL decomposes each mesh cell into sim- 
plices (e.g., each hexahedron into tetrahedra) for cell-based operations such as interpo- 
lation (Chapter 1 1). 

int res; 

float field_value, max_value = 0.0; 

FEL_vector3f coords [8]; 

FEL_c ell max_c ell; 

In addition to meshes and fields, FEL provides a number of fundamental data types, 
such as vectors (Chapter 3), positional classes (Chapter 6), and cells (Chapter 5). One 
can also declare arrays these types. For instance, the coords array above can be used 
to store the coordinates for each vertex of a cell. 

FEL_cell_iter iter; 

for (mag_grad_p_f ield->begin (&iter) ; ! iter . done (} ; ++iter) { 

FEL supports a mesh-independent way to loop over the cells of a mesh: iterators . 

Using an iterator, one can write algorithms that work with fields based on a variety of 
mesh types, without having to provide conditional statements that are type-dependent 
(Chapter 16). . 

mesh->coordinates_at_cell ( *iter , coords) ; 
centroid. set (0 . 0 , 0.0, 0.0); 

for (int i = 0; i < ( *iter) .get_n_nodes { ) ; i++) 

centroid += coords [i]; _ _ - 

centroid /= ( *iter) .get_n_nodes ( ) ; 


1.2. A GRAPHIC EXAMPLE 
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A fundamental operation in FEL is querying a mesh for geometric information, in 
this case the coordinates of the cell currently represented by the iterator. The coor- 
dinates_at_cell call returns coordinates for each vertex of the cell (Chapter 1 1). 

The cell, in turn, is queried for its number of nodes in order to calculate the centroid. 

res = n\ag_gracLp_f ield->at__phys_pos (centroid, &f ield__value) ; 
if (res ! = FEL_OK) continue; 

cout << "Pressure gradient magnitude at " « *iter 
<< " is " << field_value << endl; 
if (field_value > max_value) { 
max_ value = field_value; 
max_cell = *iter; 

} 

} 

cout << "Maximum of " << max_value << " in " « max_cell << endl; 

} 


One key group of member functions on FEL fields are the “at” calls, which support 
requests for field values. Some types of at calls merely retrieve field values at nodes, 
while others, such as at_phys.pos support field value queries at arbitrary physical 
positions. To comply with the at-phys.pos call in the example, FEL must find the 
mesh cell containing (centroid) and then use the geometry of the cell and the field 
values at its vertices to interpolate and get the field-value at centroid. FEL 
provides alternate versions of at-phys_pos (via function overloading) where the user 
can provide extra arguments intended to accelerate the point location (Chapter 1 1) or 
interpolation (Chapter 17) steps. 

— — -—The at call returns 1 to signify success. It is possible that the at call in the example 
could fail in some cases. For instance, if the cell containing centroid were at the 
boundary of the mesh and had nonplanar faces, then the point location code might 
conclude that the given point is outside the mesh. In general, it is important that the 
user check the return value of the at calls in order to be assured that the field values 
produced by the call are valid. 

This example is representative of the sorts of things one can do with FEL. The re- 
mainder of this document provides further details, options, examples, and possibilities. 
The accompanying FEL Reference Manual provides a list of the FEL classes and a 
summary of the public member functions for each class. 


1.2 A graphic example 

FEL per se does not do graphics, but some of the first applications written with FEL 
do. Graphics are fun, so we include one figure illustrating the output of a visualization 
tool called gel, which was written with FEL. Figure 1.1 is just a small taste of what 
one can do with FEL; there will be more pictures to follow. 
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Figure l .I: A space shuttle launch vehicle visualization: the left side of the model 
displays ediies from the mesh, the right side shows contour lines lor a pressure derived 
field, i Data hv Pieter Bunina visualization courtesy of Tim Sanclsfrom.) 
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1.3 FEL1 and FEL2 

An initial version of the Field Encapsulation Library was first presented at Visualiza- 
tion ’96 [BKGY96]. The library described in this document represents a fundamental 
redesign and complete rewrite of “FEL1”. Some of the features that are new to the 2.0 
release of FEL include: 

• derived fields, with demand-driven evaluation 

• built-in support for the over 50 standard derived fields defined by 
PLOT3D [WBPE92] 

• differential operators, employing either first- or second-order techniques, also 
with demand-driven evaluation 

• transformed meshes and fields 

• iterators 

• support for arbitrary field node type, based on C++ templates 

• an interface supporting the easy composition of fields 

• much better opportunities for code reuse in the development of new classes 

• support for a variety of cell types 

• multiple, user-selectable spatial interpolation modes 

• multiple, user-selectable simplicial decomposition modes 

• working set management of time-series data 

• temporal interpolation 

• significantly improved memory management, including reference counting for 
meshes, fields, and interpolants 

• support for more mesh types, including regular and unstructured 

• support for surface meshes 

• additional file reading capabilities, including support for reading PLOT3D 
“FORTRAN unformatted” files on either a workstation or Cray 

• support for the reuse of the interpolation data associated with a cell 

• integrated support for paged meshes and fields 

1.4 Acknowledgements 

We would like to thank the members of the NAS Data Analysis Group for their con- 
structive feedback, patience, and encouragement during the FEL2 development pro- 
cess. 
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Chapter 2 

Templates and Tvpedefs 


FEL makes significant use of two C++ features-templatesand typedefs-in order to in- 
crease the flexibility and power of the library without decreasing ease of use. Typedefs 
are available in C, thus they are likely to be familiar to many users. Templates, on the 
other hand, are new with C++ and are probably less familiar to most. In FEL the two 
features are used together in order to make it easy for the user to reap the benefits of 
templates without having to learn template syntax. 


2.1 Templates 

Templates are a feature of C++ allowing one to implement algorithms and classes 
where the specification of particular types is deferred. For instance, one could have 
a linked list of objects where the object type is specified by a token “T”. The user 
can then declare that he or she wants a list of objects of a particular type, for instance 
float, by using the template syntax: list<f loat>. Informally, one can imagine 
the compiler and linker working together to generate the code required by a particular 
instantiation, such as a list of floats, by outputing customized source with the global 
substitution float for T, and then compiling the dynamically generated source. How 
and when template code is actually instantiated varies from system to system. The user 
experiences these differences as variations in how programs are compiled and linked, 
but the variations do not change the style in which one writes applications using FEL. 

The most prominent use of templates in FEL is in the node type T of 
FEL-typed-f ield<T>. Templates allow the library to use one common set of 
classes for representing fields, regardless of field node type. There are two major ben- 
efits of the templated node approach. The first benefit has implications primarily for 
the user. The use of templated nodes means that one can create fields with node types 
which are not built into the library, without modifying the library. For instance, a com- 
putational scientist could create a scalar field where each scalar is represented by a dou- 
ble by writing FEL.core_f ield<double>. (The class FEL.coreJ: ield<T> is 
described in Chapter 18). The compiler and linker automatically handle generating the 
appropriate code. One could also instantiate fields where the node type is a structure 
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containing an aggregation of values, such as the solution vector used by a particular 
solver. 

A second benefit of the template approach impacts the developers of FEL more di- 
rectly than the users. The templated field node type allows the consolidation of classes, 
such as those for scalar fields and vector fields, which were formerly separate. This 
consolidation has many positive benefits: templated fields mean that there is less li- 
brary code to maintain, and bug fixes in the field classes do not have to be propagated 
to fields of each node type variation. 


2.2 Typedefs 

Typedefs (provided by both C and C++) allow one to define convenient names for what 
are typically more complicated types. In FEL, typedefs are used primarily to hide 
template syntax from the users. For example, the statement: 

typedef FEL_typed_f ield<f loat> FEL_f loat_f ield; 

makes it possible to write FEL_f loat-f ield where one would have to write 
FEL_typed_f ield< float > without the typedef. Further typedefs are provided for 
working with FEL^typed-f ielckf loat> instances. For example, the statement: 

typedef FEL_pointer<FEL_typed_f ield< float > > 

FEL_f loat_f ield_ptr ; 

makes it easy to declare a pointer to a float field (FEL_f loat.field.pt r) without 
using the more verbose template syntax. FEL provides typedef names for the most 
commonly used template instantiations; in particular, the types related to float fields 
and vector fields (where the vectors are composed of 3 floats) are included in the li- 
brary. For many applications, one can use FEL without ever using template syntax. For 
more advanced users, the option of going to template notation is always available. For 
example, to instantiate a field with a new node type would require using at least some 
template notation. 

The typedef names built into FEL are described in the following chapters, intro- 
duced with the templated classes for which they are written. 



Chapter 3 

Vector and Matrix Classes 


FEL provides basic vector and matrix support for the user The support is primarily 
for small vectors (2, 3, or 4 components) and small, square matrices. The vector and 
matrix classes are written using C++ templates, and typedefs are provided for the most 
commonly used instantiations. 


3.1 The vector classes 

The vector classes are as follows and are listed in Table 3.1. FEL provides specific 
classes for vectors of length 2, 3, or 4, rather than using the arbritrary length vector 
class FEL_vector<N, T>, for efficiency reasons. 

The vector class typedefs built-in to FEL are listed in Table 3.2. The convention 
for typedef names is to append the original templated class name by the vector length, 
if not already provided, and a letter specifying the component type. The suffixes used 
in FEL are i, f, and d, corresponding to the C types int, float, and double, 
respectively. 

FEL provides most of the basic math operators for the vector classes, so for exam- 
ple, one can write statements such as: 

float f; 

FEL_vector3f a, b; 

FEL_vector3f c(1.0, 2.0, 3.0), d<13.0, 17.3, 21.1); 


Class 

Description 

FEL-vector2<T> 

vector of 2 type T components 

FEL_vector3<T> 

vector of 3 type T components 

FEL.vector4<T> 

vector of 4 type T components 

FEL_vector<N, T> 

vector of N type T components 


Table 3.1: The basic vector class templates. 
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Typedef 

Template Class 

Component Type 

FEL_vector2i 

FEL_vector2<T> 

int 

FEL_vector2f 

FEL.vector2<T> 

float 

FEL_vector2d 

FEL.vector2<T> 

double 

FEL.vector3i 

FEL.vec t or3 <T> 

int 

FEL_vector3f 

FEL.vector3<T> 

float 

FEL_vector3d 

FEL.vector3<T> 

double 

FEL^vector4i 

FEL-vector4<T> 

int 

FEL^vector 4 f 

FEL.vector4<T> 

float 

FEL-vector4d 

FEL_vector4<T> 

double 


Table 3.2: The vector typedefs. 


a = 3 * c + d; 

b = a + c - d; 

f = a [ 0 ] ; 

b [ 0 ] = 10.0; 

b. set (10.0, 20.0 30.0); 

cout << "a = " << a << endl; 

The final statement shows the use of the C++ style output stream operator which 
can be handy for writing out the contents of vectors, for example for debugging. See 
the file FEL.vector .h for a complete listing of the operations supported for vector 
objects. 

3,1,1 FEL vectors as arguments to other libraries 

Occasionally, one may wish to use FEL vector values as arguments to routines provided 
by other libraries. For example, one may want to make calls to OpenGL graphics 
library routines taking vector arguments, without having to repackage the data into a 
different structure. The FEL vector method v ( ) returns a pointer which can be used 
with routines requiring a C-style pointer to the beginning of an array. For instance: 

FEL_vector3f v3f(1.0, 2.0, 3.0); 

FEL__vect:or3d v3d(4.0, 5.0, 6.0); 
glVertex3fv(v3f . v{ ) } ; 
glVertex3dv ( v3d . v ( ) } ; 

The v() method is defined in FEL,vector.h and should be easily inlined by 
most compilers. 


3.2 The matrix classes 

FEL provides typedefs for the instantiations for float and double versions of the N x N 
matrices, where N is 2, 3, 4, 5, 6, or 8. The suffix convention follows along the same 
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Class 

Description 

FEL-matrix22<T> 

2 by 2 matrix of type T components 

FEL_matrix33<T> 

3 by 3 matrix of type T components 

FEL_matrix44<T> 

4 by 4 matrix of type T components 

FEL-matrix55<T> 

5 by 5 matrix of type T components 

FEL_matrix66<T> 

6 by 6 matrix of type T components 

FEL_matrix88<T> 

8 by 8 matrix of type T components 


Table 3.3: The matrix templates. 


lines as that used for vectors. For instance, FEL_matrix33f signifies a 3 by 3 matrix 
with float components. Internally matrices are stored in row-major order, i.e., “C” 
style rather than FORTRAN style. FEL provides the basic math operators for matrices 
and for matrices with vectors, such as multiplication. The library also provides the 
function FEL_invert for inverting matrices. The inversion code is written using an 
analytical technique for FEL_mat:rix22<T> and FEL-matrix33<T>. For larger 
matrices the library uses Gauss-Jordan elimination with partial pivoting. See the file 
FEL_matrix . h for a complete listing of the operations supported for matrix objects. 
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Chapter 4 

The FELtime Class 


FEL represents time using the class FEL_time. A class is used to represent time, 
rather than simply a C float or int, because there is more than one time representation 
in which a suer may want to work. FEL-time supports 3 representations: 

• physical 

• computational 

• integer computational 

Physical time corresponds to the non-dimensional time used by the flow solver in an 
unsteady flow simulation. Computational time corresponds to a temporal discretiza- 
tion, where a time-varying data set is represented as a sequence of time steps. The 
first time step would be at computational time 0, the second at 1, and so on. The user 
can request data at a computational time intermediate to the time steps by providing a 
value with a fractional part; e.g., computational time 0.5 to get data temporally inter- 
polated half way between the first and second time steps. Integer computational time, 
also known as “time step”, specifies a specific step in a time-series data set, without 
temporal interpolation. 

FEL.t ime contains a tag indicating the representation currently in use and a union 
of an int and a float where the time itself is stored. By default the time representation is 
FEL -TIME-RE PRESENTATION_UNDE FINED. The user can set time using calls such 
as set-physical. See FEL-time.h for the list of set_* and get-* calls. In 
most applications, one will typically work with one representation for time. Physical 
time provides a representation that is independent of how a flow solver chose to save 
snapshots through time, e.g., whether the data were written out at a fixed physical 
time interval or at adaptive intervals. Working in terms of time steps makes it easy 
for the user to access the data in a time series without regard to how the data are 
positioned in physical time, e.g., to compute statistics over all the time steps. Floating- 
point computational time supports the added flexibility of expressing times requiring 
temporal interpolation. 
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4.1 Multiple time representations: A caveat 

In applications where one works with both physical and computational time simultane- 
ously, there is a caveat. The class FEL_time defines the operator ==, i.e., an equality 
test. When objects containing time, such as physical positions and vertex positions 
(described in the following chapters), are compared for equality, they compare their 
respective time values. Given a physical time value and a computational time value, 
an FEL-time object would have to convert one representation to the other In order 
to compare like representations. Unfortunately, FEL.time objects do not contain the 
information necessary to map between physical and computational time. In general the 
choice of which mapping to use is ambiguous, since there could be many time-varying 
fields, and they may not all use the same mapping. FEL-time operator == returns 
false if it cannot compare the times. 

Mixing time representations in the arguments to FEL calls should not cause prob- 
lems internal to FEL, to the best of our knowledge. Mixing time representions when 
doing at.phys.pos calls may cause some minor inefficiencies, though it should 
not cause the library to return incorrect answers. This case is revisited later when 
we describe the use of at_phys_pos (Chapter 17). The time mixing caveat should 
only be of concern to the user if he or she uses the == operator with the positional 
classes described in the next chapter: FEL_vertex-cell, FEL.phys.pos and 
FEL_s true Cur ed.pos. 

There is no easy solution to the “got one representation, need another, don’t know 
how to convert” problem. The library could choose one time format to use throughout 
the application programmer’s interface, but no matter what the choice, there will be 
cases where the selected format is awkward for the user. Physical time is a higher-level 
representation, but it is not convenient when one wants to loop over time steps. On the 
other hand, some applications would like to treat an unsteady data set as continuous 
in time and are not interested in how the samples were chosen, thus physical time is a 
better choice than computational (or time steps). The FEL design opts for flexibility 
and ease-of-use, at the risk of unexpected behavior in some obscure cases. 


Chapter 5 


The FEL_cell Class 


The FEL-cell class defines a general purpose object for representing cells, the basic 
discretization building block for computational domains. FEL uses a general cell def- 
inition, where cell types include vertices, edges, triangles, quadrilaterals, tetrahedra, 
and hexahedra. Cells are used, for example, for point location: given a physical point 
and a hexahedral mesh, the library can return a cell containing the point. The cell types 
which can be represented by FEL-cell are listed in Table 5.1. 

Cells may be grouped by their dimension and referred to as fc-cells or k- 
dimensional cells. Every cell type has an associated dimension, the Dimension column 
of Table 5.1 lists the dimension for each cell type supported by FEL. Writing in terms 
of A;-cells is more convenient in some cases, for instance, one can speak of the pri- 
mary cells of a surface being 2-cells, without having to distinguish between triangles, 
quadrilaterals, or some other type of polygon. 

An important concept related to cells is faces. A cell c in a mesh M. is the face 
of another cell d in M if c is defined by a subset of the vertices defining d. The most 
common use of the term face is with hexahedral and tetrahedral cells: a hexahedron has 
quadrilateral faces, and a tetrahedron has triangle faces. But the term is more general: 
hexahedra and tetrahedra also have edge and vertex faces, even an edge has vertex 
faces. The topic of ceils, faces, and incidence relationships will be revisited in the 
mesh chapter (Chapter 1 1). 


Cell 'type 

Letter 

Dimension 

Vertices 

FEL .CELL .VERTEX 

V 

0 

HHHI 

FEL.CELL.EDGE 

E 

1 


FEL -CELL -TRIANGLE 

F 

2 


FEL-C ELL-QUADRILATERAL 

Q 

2 

HEH 

FEL-CELL.TETRAHEDRON 

T 

3 

MDH 

FEL-C ELL-HEXAHEDRON 

H 

3 

8 I 


Table 5.1 : The FEL cell types. 
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Name 

Type 

Default 

type 

FEL.cell-type.enuin 

FEL.CELL.UNDEFINED 

subid 

short 

-1 

ijk 

FEL.vector3i 

(none) 

zone 

int 

FEL.ZONE.UNDEFINED 

time 

FEL.time 

FEL_time ( ) 


Table 5.2: The FEL.cell data members. 


5.1 The cell data members 

The 5 data members of an FEL.cell object are listed in Table 5.2. The type speci- 
fies one of the types listed in Table 5.1, or FEL.CELL.UNDEFINED. The zone would 
specify a zone in a multi-zone mesh. Time is used when working with time-varying 
data. The ijk and subid data members are used to specify a particular cell of a given 
type. The particular usage of ijk and subid for each type of cell depends on whether 
the mesh is structured or unstructured. See Chapters 12 and 13. 


5.2 The FEL_vertex_cell class 

An FEL.ver tex.cel 1 object represents a vertex in a mesh. The class is derived from 
FEL.cell and can be used as an argument to any routine expecting an FEL.cell ar- 
gument. The FEL.vertex.cell class exists in order to make the development of 
algorithms written in terms of mesh vertices a bit easier. The class provides the oppor- 
tunity to get some compile-time checking of routines which work only with vertices; 
and FEL provides routines for vertex cells which have some slight optimizations over 
the general cell routines. 




Chapter 6 

Positional Classes 


The position of an individual point in space can be represented by an 
FEL-ver t ex.cel 1, FEL_phys.pos,or FEL_structured_pos object, see Ta- 
ble 6. i . All three classes also contain a representation for time, so the objects in general 
can represent a position in space and time. 

6.1 FEL_vertex_cell 

The FEL-vertex.cell class is used for specifying an individual vertex in a mesh. 
FEL_ver t ex.cel 1 is derived from the class FEL.cell and inherits most of its func- 
tionality from the cell class. FEL.cell and FEL.vertex-cell are described in the 
previous chapter. 

6.2 FEL_phys_pos 

An FEL.phys.pos object represents a position in physical space and time. One can 
request field values at a physical position using statements such as: 

int res ; 
float f; 

FEL_phys_pos p(1.2f, 3.3f, O.Of); 
res = f loat_f ield->at_phys_pos (p, &f ) ; 
if (res == 1) 

cout << 11 the field at"" << p << " is " << f « endl ; 


Class 

Derived From 

FEL.ver t ex-cel 1 

FEL.ceil 

FEL-phys.pos 

FEL-vector3f 

FEL.structured-pos 

FEL_vector3f 


Table 6.1: The FEL positional classes and their parent classes. 
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Since FEL.phys.pos is derived from FEL.vector3f, FEL.phys.pos can be 
used as with routines expecting an FEL_vector3f argument; the time component 
of the FEL.phys.pos is ignored. For example, to get the physical coordinates at a 
vertex, one could write: 

mesh->coordinates_at_vertex_cell (vertex_cell , &phys_pos) 

One can also use FEL.phys.pos objects with the mathematical operators += and 
-=, where the right-hand-side argument is of type FEL.vector3f . For instance, one 
could add to the spatial coordinates of a physical position by writing: 

FEL_phys_pos p(10.0f, 11. Of, 12. Of); 
cout << "before: ” << p « endl; 
p += FEL_vector3f (1 . If , 2.2f, 3.3f); 
cout « "after: " « p « endl; 

The statements with cout illustrate the use of the C++ ostream operator defined for 
FEL.phys.pos. The ostream operator displays the physical coordinates; the time 
is also displayed if it is defined. The output from the statements above would look like: 

before: TlO, 11, 12) 
after: (11.1, 13.2, 15.3) 

6.3 FEL-struetured-pos 

The class FEL_structured-pos represents positions in structured meshes de- 
scribed by floating-point i, j, and k values. FEL-structured_pos val- 
ues are sometimes known as “computational coordinates”, though they are only 
applicable to structured meshes (or a structured zone in a multi-zone mesh). 
FEL.structured.pos objects also contain an integer zone value and time. The 
zone comes into play when working with FEL_multi_inesh objects, and time is rel- 
evant to time- varying objects. 


6.4 FEL_vector3f_and _int 

The FEL_vector3 f .and-int class, as its name suggests, represents objects con- 
sisting of a vector of 3 floats and an integer. The FEL_vector3 f _and_int 
class is used primarily for working with meshes that include what is known 
in PLOT3D data sets as IBLANK [WBPE92]. IBLANK values are inte- 
gers, one per vertex. Users see the FEL-vector3f-and_int type in the 
mesh calls coordinates.and-iblank_at-vertex-cell ( ) and coordi- 
nates_and-iblank_at_cell { ) . See Chapter 1 1 for more on these routines. 


Chapter 7 

Memory Management 


FEL allocates and deallocates memory using the standard C++ new and delete oper- 
ators. In general, FEL objects that allocate memory are also responsible for deallocat- 
ing the same memory. Likewise, memory allocated by the user is the user’s responsibil- 
ity for deallocation. There is one key exception to this convention with FEL: memory 
allocated for buffers passed to constructors of mesh and field subclasses becomes the 
responsibility of FEL to deallocate. For instance, a buffer with float values passed to 
the constructor of an FEL-core-f loat.f ield would become the responsibility of 
the core field to deallocate when the core field is destructed. Users who work solely 
with meshes and fields constructed via the file reader functions (see Chapter 22) do not 
have to be concerned with this issue, since the the convenience functions allocate and 
pass the buffers appropriately. If the user constructs a mesh or field explicitly, providing 
data buffer arguments, then it necessary that the buffers be allocated using the C++ new 
[ ] operator, i.e. the C++ operator for allocating arrays of objects. This requirement 
is necessary because ultimately FEL will deallocate the buffer using the C++ delete 
[ ] operator, and in general it is hazardous to" allocate memory in one manner (e.g. 
malloc) and deallocate in another (e.g. delete [ ] ). With FEL.core.f ield in- 
stances, the user also has the option of providing a flag that will suppress the field’s 
attempt to deallocate the solution data buffer. In this case, how the buffer was allocated 
is irrelevant to FEL, and the ultimate responsibility for deallocating the buffer would 
remain with the user. 


7.1 Reference counted objects and pointers 

Reference counting is a technique for tracking the number of references to each object 
in a set of objects. Reference counting supports sharing objects and the detection of 
a safe time to deallocate an object (i.e., when the number of references to an object 
goes to 0). In FEL, the most prominent use of reference counting is for mesh and field 
instances. Reference counting allows the user to create meshes and fields and then 
build additional fields, such as derived fields and differential operator fields, in terms 
of the original fields. The reference count book-keeping is done transparently. The 
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reference counting frees the user from some relatively tedious and error prone work 
and in general is unobtrusive. 

FEL uses two classes that work together to implement reference counting: 
FEL_ref erence_counted_obj ect and FEL_pointer<T>. The FEL_mesh 
and FEL_f ield classes inherit from FEL_reference_counted-object. The 
FEL-pointer<T> class is used to represent references to an object derived from 
FEL_ref erence.counted.obj ect. The template type T stands for the particular 
class being pointed to, e.g„ FEL_mesh. FEL provides typedef names for most point- 
ers to reference counted objects, using the convention that pointer names take the class 
name of the object pointed to and append For example, FEL_mesh_ptr and 

FEL-f loat.f ield-ptr represent pointers to meshes and float fields. The pointer 
class provides definitions for the operations that can potentially change the reference 
count of an object. For instance, FEL_pointer<T> defines an assignment operator, 
so that when one pointer is assigned to another, one reference is potentially incremented 
and one is potentially decremented. We say "potentially” since pointers can have the 
value NULL, and no counts need to be changed when a pointer has a NULL value. 

The FEL_pointer<T> class is designed so that one can use pointer instances in 
a manner similar to raw C-style pointers, e.g., one can write statements using -> 
syntax. There are a few cases where FEL pointer instances differ in usage from raw 
pointers: 

• default initialization 

• casting 

• comparisons with NULL 

• deleting 

FEL pointers, unlike raw C-style pointers, are guaranteed to always be initialized to 
NULL by default. This difference is not necessarily noticeable to the user, since one 
can still manually initialize pointers as one would do with raw pointers. 

A second case where FEL-pointer<T> instances differ from raw pointers is in 
downcasting. For example, one may have a pointer to an FEL_mesh, but what one may 
need is a pointer to an FEL.structured.mesh, in order to access a method that is 
specific to structured meshes. With raw pointers one could write: 

FEL_mesh* mesh; 

FEL_s t rue tur ed_mesh* s t rue tured_mesh ; 

structured_mesh = (FEL_structured_mesh*) mesh; 

Unfortunately, if FEL_mesh* and FEL.structured_mesh* were replaced by 
FEL.mesh.ptr and FEL.structuredjmesh.ptr, the previous excerpt would not 
work. The cast fails because FEL.mesh.ptr (i.e., FEL.pointer<FELjmesh>) 
technically is not a C pointer; there is no direct conversion from FELunesh.ptr to 
FEL-Structured-mesh-ptr (i.e., FEL.pointer<FEL.structuredjnesh>). 
The cast can be accomplished, but the required syntax is awkward. For the previous 
example: 
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FEL_mesh_ptr mesh; 

FEL_structured_mesh _ptr structured_mesh; 

structured_mesh - (FEL_structured_mesh*) (FEL_mesh*) mesh; 

The less-than-obvious syntax is one motivation for providing downcasting macros for 
the user. See the following chapter for a list of the casting macros provided by FEL 
and some additional motivation for using the macros. 

A third case where FEL_pointer<T> instances differ from raw pointers occurs 
when testing whether or not a pointer has the value NULL. FEL_pointer<T> class 
has the member function null ( ) that can be used for this purpose, for example: 

FEL_mesh_ptr mesh; // default initialization 

assert (mesh .null ( ) ) ; 

. „ . assign non-NULL value to mesh 

assert ( Imesh . null ( ) ) ; 

Two FEL.pointer<T> instances (both instantiated with the same type for T) can 
also be tested for equality or lack there of using the == or ! = operators. 

FEL pointer objects differ from raw pointers in a fourth respect: deletion. Whereas 
in genera] one can call delete directly on a heap allocated object, the equivalent 
statement will not work with an FEL pointer instance. In a sense, this should not 
be surprising, reference counting takes responsibility for deallocating an object when 
there are no more references; the user should not have to take the same responsibility. 
Nevertheless, there are a few occasions where the user may want to cause the deallo- 
cation to happen sooner. For example, one may want to free memory used by a large 
mesh or field that one knows will not be used again. One can do this indirectly by 
assigning NULL to FEL pointers. Assuming that there are no other references to the 
mesh or field, assigning NULL to the remaining reference will decrement the count to 
0 and induce the reference counting mechanism to do the deallocation. 

Finally, some readers might observe in the previous discussion that it is possible 
to obtain the raw pointers to reference counted objects. One may be tempted to do 
this, since users are more familiar with C-style pointers and casting than the refer- 
ence counting support classes of FEL. One should, however, resist this temptation. 
Reference counting requires little overhead, so performance is typically not an issue. 
Furthermore, the reference counting mechanism provided by FEL is robust, but all bets 
are off if one breaks through the abstraction and starts handling raw pointers directly. 
The potential consequences for defeating the reference counting are memory misman- 
agement bugs which can be very difficult to track down. 


7.2 Reference counting and mutual exclusion 

The FEL reference counting mechanism is also designed to support the use of reference 
counted objects in a multi-threaded environment. In a multi-threaded scenario, it is pos- 
sible that several threads can take actions that change the reference count of a shared 
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object, such as mesh or field. Such actions include assigning one FEL_pointer<T> 
instance to another or passing a pointer object to a subroutine. Since such actions can 
take place simultaneously, it is important that the changes to the reference counts be 
made atomically. In other words, only one thread should be able to change the count 
at one time. FEL provides reference counting with critical section protection for count 
changes via the class FEL_mutex_ref erence-COuntecLobj ec t which is derived 
from FEL-ref erence-counted_ob j ec t. The ’’mutex” version of the class uses 
a locking primitive provided by the task synchronization library (TSLonutex) to en- 
sure mutual exclusion. Meshes, fields, and interpolants in FEL are all derived from 
FEL_mutex_ref erence.countecLobj ect, thus the reference counting for those 
objects should still work reliably in multi-threaded applications. 


Chapter 8 

Dynamic Casting 


There are occasions when one has a pointer to a base class, but what one needs is a 
pointer to a class derived from the base class. For example, one may have a pointer 
to a mesh: FEL-mesh_pt:r, but need a pointer to a structured mesh (a subclass of 
mesh): FEL.s true turedjnesh.pt r. The typical case where a pointer to a subclass 
is necessary arises when one needs to call a method that is available only from the 
subclass, for instance a method specific to structured meshes. Casting from one type to 
another type that is lower in the same class hierarchy is known as downcasting . FEL 
is designed to try to minimize the amount of downcasting required, nevertheless there 
are still cases where it is necessary. One hazard of downcasting is that it is possible 
to do it incorrectly, e.g., if the mesh pointer were referring to an unstructured mesh, a 
structured mesh cast would lead to dire consequences. Unfortunately, incorrect casts 
in general cannot be detected at compile-time; it is not until run-time that they become 
apparent. 

C++ provides a standard way to downcast safely: dynamic casts. Using a dynamic 
cast, one can downcast a pointer to a particular subclass. If the cast is legal, then the 
result is a pointer to the subclass; if the cast is not legal, then the result is NULL. Thus 
one can downcast and then check if the cast completed successfully before continuing. 
Dynamic casting is a relatively new C++ feature, and it is not supported by some older 
compilers. FEL provides macros that behave like dynamic casts for the most com- 
mon casts of FEL objects. The naming convention for the casts is FEL/_TO_f-CAST, 
where an FEL pointer to a type / cast to an FEL pointer to a type t. For example 
FEL-MESH.TO.STRUCTURED_MESH.CAST { ) takes an FEL_mesh.pt:r argument 
and returns an FEL.structured-xnesh.ptr, or NULL if the cast is not legal. One 
difference between the FEL macros and C++ dynamic casting comes in the case where 
the argument is a NULL pointer. The C++ standard dictates that the dynamic cast 
of a NULL pointer should cause an exception, the FEL macro simply returns NULL. 
Table 8.1 lists the casting macros provided by FEL. Though there are many macros 
defined by the library, the need for them in user code is relatively infrequent. As the 
design of the library evolves, we are working towards further reducing the need for 
downcasting calls. 
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Cast 

FEL-OBJECT-TO-MESH.CAST ( ) 

FEL-OBJECT.TO-FIELD.CAST ( ) 

FEL J^ESH.TO -STRUCTURED J4ESH.CAST ( ) 

FEL-MESH-TO.UNSTRUCTURED -MESH-CAST ( ) 

FEL_F I ELD-TO.FLOAT-F I ELD-CAST ( ) 

FEL-F I ELD-TO.VECTOR3 F-F I ELD-CAST ( ) 

FEL-FIELD-T0-PL0T3D-Q.FIELD.CAST ( ) 

FEL-FIELD-TO-CORE-FLOAT -FIELD-CAST ( ) 

FEL-FIELD-TO-CORE.VECTOR3F-FIELD-CAST ( ) 

FEL-FIELD-TO-TIME-SERIES-FLOAT .FIELD-CAST ( ) 

FEL-FIELD-TO-TIME.SERIES-VECTOR3F-FIELD-CAST ( ) 

FEL-FIELD-TO-TIME-SERIES-PLOT3D-Q-FIELD.CAST ( ) 

FEL-INTERPOLANT-TO-HEXAHEDRAL-I SOPARAMETRIC -CAST ( ) 


Table 8.1: The FEL dynamic casting macros. 


User code with an FEL dynamic cast macro would typically look like: 

FEL_mesh_ptr mesh; 

FEL_structured_mesh_ptr structured_mesh; 

structured_mesh = FEL_MESH_TO_STRUCTURED_MESH_CAST (mesh) 

if ( ! s t rue tured_mesh. null { ) ) { 

} 

else { 

error . . . 

} 

There is a second reason for the use of macros with downcasting. Pointers to FEL 
meshes and fields are not raw C-style pointers; rather the types such as FEL_mesh_ptr 
or FEL-f loat.f ield-ptr are typedef names for FEL-pointer objects. The 
classes FEL.pointer and FEL_ref erence.counted.object are used together 
to provide reference counting support in FEL (as described in the previous chapter). 
Unfortunately, straight-forward C-style downcasting does not work with FEL pointers, 
though it is still possible to downcast using a more arcane syntax. The macros are 
intended to hide that syntax. 



Chapter 9 

Interpolation 


One of the most useful features of FEL is that it allows gridded, spatially discrete data 
to be treated as a continuous domain. In fact, this capability supports the key abstrac- 
tion of a field and is one of the primary aims of the library. The central mechanism 
underlying the field abstraction is spatial interpolation. FEL supports temporal inter- 
polation, too; this is described separately in Chapter 25. 

9.1 Setting interpolation modes 

In FEL, queries for field values at arbitrary physical space locations (i.e., 
at-phys.pos ( ) , see Chapter 17) invoke a point location algorithm which finds the 
mesh cell containing the query point. FEL then uses the geometry of this cell, and the 
field values at its vertices, to determine a field value at an interior point. The last step 
can be done in one of three ways, and you tell FEL which method you want with a set 
command. 

FEL_mesh_ptr mesh; // defined elsewhere. . . 

mesh->set ( FEL_INTERPOLATION, FEL_NEAREST_NEIGHBOR_INTERPOLATION) 
mesh->set (FEL_INTERPOLATION, FEL_ISOPARAMETRIC_INTERPOLATION) ; 
mesh->set (FEL_INTERPOLATION, FEL_PHYSICAL_SPACE_INTERPOLATION) ; 

These set calls are really methods on FELjnesh, but you can also call them on 
any field, and it will simply forward the call to its mesh. Note that however it is set, 
the interpolation mode on a mesh affects all fields based on that mesh. This rather 
unfortunate state of affairs will be amended in a future release of FEL so that fields will 
be less dependent entities. In the meantime, if you want different interpolation modes 
on fields sharing a common mesh, set the interpolation mode just prior to querying 
the field — but be aware that this may not be a robust strategy in a multithreaded 
application. 

The interpolation modes work in concert with the simplicial decomposition mode, 
since the cell type which the interpolation routines receive may be altered by simplicial 
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decomposition. (Details on simplicial decomposition may be found in Chapter 12.) 
Thus the three interpolation modes listed above, combined with the three choices for 
simplicial decomposition — 

// setting decomposition off, even and odd: 
mesh->set (FEL_SIMPLICIAL_DECOMPOSITION, 0) ; 
mesh->set (FEL_SIMPLICIAL_DECOMPOSITION, 1) ; 
mesh->set (FEL_SIMPLICIAL_DECOMPOSITION, 2) ; 

— form a 3 x 3 matrix of possibilities. If you don’t explicitly set them, the interpolation 
mode defaults to FELLESOPARAMETRIC-INTERPOLATION, and “even” simplicial 
decomposition is turned on (FEL-SIMPLICIALJ3ECOMPOSITIQN, 1). The dif- 
ferent combinations of interpolation and simplicial decomposition may produce quite 
different numerical results on the same dataset, and they vary widely in the amount of 
numerical work involved. A few general considerations about choosing these modes 
are presented below (and see [KL95] for a comparison of isoparametric and physical 
space interpolation methods on tetrahedra), but for the most part, the choice is context 
dependent. 

The combination of interpolation and decomposition modes affects the differential 
operator fields, since they obtain their required spatial derivatives by analytically dif- 
ferentiating the interpolating polynomials. This is the case even when the differential 
operator fields are queried at a vertex. For more information on the differential operator 
fields, see Chapter 20. 


9.2 Nearest neighbor interpolation 

FEL-NEARESTJJEIGHBOR-INTERPOLATION really doesn’t do interpolation, 
rather it assigns to an interior point the field value of the nearest vertex of the en- 
closing cell. This “interpolation” method is the fastest of the bunch, but obviously not 
very smooth. 


9.3 Isoparametric interpolation 

FEL-ISOPARAMETRIC -INTERPOLATION transforms the enclosing cell and the 
query point into computational space and then interpolates a field value at the query 
location using linear basis functions. If the cell is a tetrahedron, the physical to compu- 
tational space transformation of the query location is done analytically. For hexahedra, 
however, the physical to computational space transformation of die query location must 
be done numerically, using Newton’s method, and this requires evaluating and inverting 
the Jacobian matrix for each iteration. 

In computational space, the coordinates of the query point are f , r C* and the field 
value / = /(£, 77 , Q. Field values at cell vertices are referred to as /„, where n ranges 
from 0 to 1 less than the number of nodes in the cell. 

For tetrahedra, we use the linear basis function: 

/(£> 0 = (i*o - C “ v - O/o + f A + *7/2 + C/3 
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and for hexahedra: 
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9.4 Physical space interpolation 

FEL.PHYS ICAL.S PACE-INTERPOLATION solves for the coefficients of an interpo- 
lating polynomial in physical space, using the locations and field values of the enclos- 
ing cell vertices as constraints. The calculation does not involve any physical to com- 
putational space transformations, but it does require inversion of a 4 x 4 (tetrahedron) 
or 8 x 8 (hexahedron) matrix. 

In physical space, the coordinates of the query point are x, y, z, and the field value 
/ = f(x, y, z ). Coordinates of cell vertices are referred to as i„, y n , z n , and field val- 
ues at the corresponding vertices are referred to as / n , where n ranges from 0 to 1 less 
than the number of nodes in the cell. The coefficients of the interpolating polynomial 
are referred to as c n . 

For tetrahedra, we use the interpolating polynomial: 

/(x,y,z) = co + cix + c 2 y 4- c 3 z 
and determine the coefficients by: 
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For hexahedra, we use the interpolating polynomial: 

/ (x, y,z) = eg + cix + c 2 y + c 3 z + e 4 xy 4- c 5 xz + c 6 yz + c 7 xyz 
and determine the coefficients by: 
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The 4x4 and 8x8 matrices are inverted numerically, using Gauss- Jordan elim- 
ination with partial pivoting. We then check the “standard” matrix condition number 
using the infinity norm (|| A||oo ||-4“ J ||oo). and if this exceeds le 6 we perform a singular 
value decomposition, zero the singular values < le“ 3 and generate a pseudoinverse. 
We found the SVD to be necessary, especially in the 8 x 8 case, as many common 
data sets yield ill-conditioned matrices surprisingly often. The SVD routines will be 
made more flexible in a future release of FEL, with user-settable options replacing the 
hardcoded invocation and parameters. See, for example, [PTVF92] to learn more about 
SVD and related numerical issues. 

Although the matrix inversions involved in the physical space interpolation routines 
are expensive, especially if the SVD routines are invoked, the inverted matrices are 
determined purely by the geometry of the cell. They can be reused for successive local 
interpolations on the same or a different field if the successive queries all fall in the 
same cell. The inverted matrix is cached in the FEL.cell-interpolant which can 
be resubmitted as part of an at-phys.pos ( ) query. See Chapter 17 for details on the 
several variants of-at_phys_pos { ) . 



Chapter 10 

Meshes and Fields 


The heart of FEL is composed of mesh and field objects; the majority of the following 
chapters of this document will be on topics related to these two key types. Most of the 
common interface for meshes is defined by the class FEL_mesh (Chapter II). Most of 
the interface for fields is defined by the classes FEL_f ield and FEL.typed.f ield 
(Chapter 17). 

One goal in the design of FEL is that applications written in terms of the standard 
interfaces should work with a variety of mesh and field types. For example, a visual- 
ization technique written for scalar fields should work just as easily with a field where 
values are computed on demand, such as a derived pressure field (Chapter 19), as with 
a field where the data are precomputed and stored in memory (Chapter 18). The same 
code should also work with many other types of fields, such as those where data are 
paged in from disk on demand (Chapter 23) or fields that vary with time (Chapter 25). 

While it is not hard to see the virtues of code reuse, it is not as easy in practice 
to design interfaces that make such reuse straightforward for the user. Even with the 
best of interfaces, it still may take some effort on the user’s part to think in more 
general terms. Not every mesh is structured, not every field is steady. In general, the 
development of truly mesh and field type- independent algorithms requires effort on the 
part of both the FEL design team and the user. Making mesh and field independent 
interfaces and algorithms a reality continues to be a learning process for us all. 


10.1 Member function style 

The member functions of the mesh and field classes follow a general style where input 
arguments come first, followed by arguments pointing to the location where results 
should be written. The input arguments typically are C++ const references, so that 
their intended use should be clearer to the user, and so the compiler may have more 
opportunities to optimize. Most member functions return an integer indicating whether 
the call was successful. We look at the return values in a bit more detail next. 
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10.2 Return values 

Most mesh and field member functions return an integer status value. The value can 
indicate success (FEL_OK) or it can indicate one of a variety of errors. The complete list 
of return values can be found in FEL-returns . h. One should use the names defined 
in FEL .returns . h rather than integer values explicitly, so that in the unlikely case 
that return codes get renumbered, one’s code will still work. One return value where 
even the best of FEL programmers have lapsed into using the integer value directly is 
FEL.OK. FEL.OK is equal to I . The interchangeability of 1 and FEL .OK is so ingrained 
in the FEL programming style that it is safe to assume that the number and symbol will 
be bound together for all time. 

We conclude this chapter by strongly encouraging the users of FEL to check return 
values. There are two main reasons why one should get into this habit. First, if one is 
serious about developing algorithms that are mesh and field type-independent, then it 
is hazardous to assume when using some member function that “this call cannot fail”. 
Even though there are already a large number of mesh and field types, it is easy to fall 
into the trap of thinking in terms of just one particular type. Even if one does consider 
all the ways that a call can fail, based on the types available in FEL today, one is still not 
completely safe. In the future there will surely be more mesh and field types introduced 
into the library. Algorithms written with careful error checking now should at least be 
able to gracefully indicate that they cannot work with a new type in the future. 

A second reason to be conscientious about return value checking is that not doing 
the checking can lead to potentially insidious bugs. Since most member functions work 
by writing their results into a location passed into the call, one always has something in 
result location, whether or not the call succeeded. In some cases it may be obvious that 
the result location contains junk, but at other times the contents may seem plausible. 
For example, the result may contain a value from a previous, successful call. Checking 
return values is the only reliable way to guard against this type of problem. 


Chapter 11 

Meshes 


Meshes are one of the two key types of objects in FEL, the other being fields. Meshes 
represent discretization of a domain into a set of cells. The cell types are drawn from 
the types represented by the class FEL.cell: vertices, edges, triangles, quadrilaterals, 
tetrahedra, and hexahedra. In an FEL mesh, it is assumed that a mesh containing a cell 
c also includes all the faces of c; so, for example, a mesh with a hexahedron c would 
also have all the quadrilateral, edge, and vertex faces of c. Meshes contain both geo- 
metric and topological information. Geometric information includes the coordinates of 
vertices and the volume of cells. Topological information includes neighbor relation- 
ships among the cells and other data about how the cells are organized. For example, a 
mesh containing a tetrahedron c can return the triangle faces or vertices of c. 

In FEL, meshes are essential to the construction of fields, since every field has a 
mesh. For fields, meshes specify the location of nodes. Nodes are the points in the do- 
main where solution values are generated by the solver or acquired by experiment. FEL 
supports vertex-centered fields, in other words fields where a node is associated with 
each vertex. There are other organizations for nodes; for instance, a hexahedral mesh 
may be “cell-centered”, i.e., a node is associated with the interior of each hexahedron, 
but such configurations are currently unsupported by FEL. 

One key responsibility of meshes is point location. Given a point p and a mesh 
containing some type of 3-cells, such as tetrahedra or hexahedra, the result of point 
location will be an integer return code and, if the location effort is successful, a cell 
containing p . The concept of point location in FEL has been generalized to meshes 
which do not contain 3-cells. For instance, FEL can represent surfaces in R 3 consisting 
of triangles or quadrilaterals (2-cells). Point location with a surface mesh returns a 
2-cell from the mesh. See Chapter 12 for details of how point location is defined 
for surfaces. Efficient point location is one of the keys in FEL to providing good 
performance overall. 

A second key responsibility of meshes is to assist with interpolation. Given the cell 
c resulting from point location, a mesh can construct an “interpolant” for use later in 
interpolation. The interpolant contains information based on the geometry of c. The 
specific information contained in an interpolant depends upon the type of interpolation 
that the user has selected. For example, if c is a tetrahedron, and the prevailing inter- 
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polation mode is isoparametric, then the interpolant would contain the basis functions 
required to linearly interpolate within c. See Chapter 9 for more detail on interpolants. 


11.1 The mesh class hierarchy 

Figure 11.1 illustrates the class hierarchy for FEL meshes. FEL_mesh inherits from the 
class FELjnutex_reference.counted_object, therefore meshes are reference 
counted in FEL, and there is critical section protection for the reference counting so that 
they can be used in multi-threaded applications. FEL_mesh specifies the interface that 
all mesh classes inherit. The mesh class also provides implementation of routines that 
are used by many of the mesh subclasses. Key subclasses of FEL_mesh are described 
in the following chapters. 

11.2 Setting and getting mesh properties 

FEL meshes support a general “set” interface for setting mesh properties: 

void set (const FEL_set_keyword_enum k f int v) ; 

The set call takes a keyword k and a value v; the complete list of keywords can be 
found in the file FEL.set-keywords .h. There is also a general “get” member 
function: 

int get{const FEL_get_keyword_enum k # int* nv, int v[], 
int zone = FEL_ZONE_UNDEFINED) const; 

The get call takes a keyword k and fills in *nv integer values into the array v. Note 
that some properties, such as the dimensions of a structured mesh, require more than 
one integer to describe. In most cases the user will know how many integers are needed 
to describe a particular property, i.e., how many values will be written into v. In those 
cases the user can pass the value NULL for the nv argument. The get function takes an 
optional final argument specifying a zone. Using the zone argument with get allows 
one to make queries of a particular zone in a multi-zone mesh. The get call returns 1 
for success or an error value otherwise. 


11.3 Simplicial decomposition 

Some algorithms work in terms of the cells of a mesh, but require that the cells be sim- 
plices, i.e., that the cells be vertices, edges, triangles, or tetrahedra. FEL supports emu- 
lating the decomposition of a mesh into simplices through what is known as “simplicial 
decomposition”. We say “emulating” because internally FEL does not change the rep- 
resentation of the mesh when simplicial decomposition is requested; it only changes 
the type of cells returned by methods such as point location. Currently only structured 
meshes contain non-simplicial cells (quadrilaterals and hexahedra), thus the following 
discussion only applies to structured mesh types. In the future other mesh types may 
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Figure 11.1: The FEL mesh class hierarchy. 
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also contain non-simplicial cells, but the interface for controlling the decomposition 
will be the same. 

In general there is more than one way to subdivide a non-simplicial cell into sim- 
plices, even if one is not allowed to introduce new vertices as part of the subdivision. 
A quadrilateral can be broken into triangles in one of 2 ways, depending upon how 
the diagonal is chosen. There are 2 ways to decompose a hexahedron into 5 tetrahedra 
and many more decompositions consisting of 6 tetrahedra. A typical requirement is 
that the cells be decomposed consistently; in other words, if there is a non-simplicial 
face shared by two cells, then the decomposition chosen for each cell must result in the 
same decomposition for the shared face. The issue of consistent decompositions in FEL 
arises with structured hexahedral meshes, since adjacent hexahedra share quadrilateral 
faces. FEL supports 2 consistent simplicial decompositions for structured meshes, to 
be described in the next chapter. FEL also provides the call 1 : 

int decomposition_cells (const FEL_cell c&, 

int* nsc , FEL_cell sc [ ] ) const; 

which takes a non-simplicial cell c as an argument and returns nsc simplicial cells sc 
that c would be subdivided into, using the prevailing simplicial decomposition setting. 
(The user can query about the prevailing decomposition mode using the get member 
function described above.) The result of decomposition-cells is undefined if 
simplicial decomposition is not on when the call is made. 

To set a particular value for simplicial decomposition, one can use the call 
set (FEL-SIMPLICIALJDECOMPOSITION, i) on a mesh, where i would have 
the value 0, 1, or 2. The value 0 signifies that simplicial decomposition should be 
turned off. Values 1 and 2 each set one of 2 alternate decompositions for structured 
meshes. 


11.4 Point location and interpolation 

FEL meshes all inherit the following interface for point location: 

int locate_close_vertex_cell (const FEL_phys_pos& p, 

FEL_vertex_cell* v) const; 

int locate (const FEL_phys_pos& p, FEL_cell* c) const; 
int locate (const FEL phys_pos& p, FEL_cell& start_cell, 
FEL_cell* c) const; 

The locate.close_vertex_cell returns a vertex v that is close to p. The vertex 
is not guaranteed to be the closest, but sometimes close is good enough. The locate 
member functions are the workhorse routines for point location. Given a point p, lo- 
cate produces a cell c. For most meshes, c is a hexahedron or tetrahedron containing 

^For those less familiar with C++ notation, the const keyword may be new. When the keyword is part 
of an argument declaration, then const indicates that the function will not modify the argument. When 
const follows the closing parentheses of a class member function declaration, then const specifies that 
calling the function will not change the state of the object. 
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p. For other types of meshes, such as surface meshes, c may be a quadrilateral cell or 
triangle cell close to p. 

The second locate routine takes an extra start-cell argument. 

For several key mesh types, such as FEL,curvi linear jnesh and 
FEL-unstructured_mesh, a point is located by “walking” from cell to cell, 
until a result cell is found. If a start.cell argument is provided, then it is used 
to initialize the walking point location routine. Providing a start cell can significantly 
improve the point location performance if the cell is close to the given point. 

An FEL,interpolant contains information based on the geometry of a cell 
used for interpolation. An interpolant is specific to a particular cell, and an 
FEL_cell,interpolant pairs a cell together with its interpolant. Since meshes 
contain cell geometric data, meshes are responsible for initializing interpolants. Ini- 
tialization is done through the member functions: 

int set_interpolant (F£L_cell_interpolant*) const ; 
int set_interpolant (const FEL_cell_interpolant & pci, 

FEL_cell_interpolant* ci) const; 
int locate_and_set_interpolant ( const FEL_phys_pos&, 

FEL_cel l_int erpolant * ) cons t ; 
int locate_and_set_interpolantTconst FEL_phys_pos&, 

FEL_cell_interpolant&, 
FEL_cell_interpolant* ) const ; 

The second set, int erpolant call provides the opportunity to reuse the interpolant 
loaded by a previous set,interpolant call, if the cell in ci is the same as in 
pci, and the interpolation mode (e.g., isoparametric) is the same. The latter two calls 
combine point location and setting an interpolant into one call. 


11.5 Coordinates 

The user can access the coordinates v of the vertices of a cell c via the calls: 

int coordinates_at_cell (const FEL_cell& c, 

FEL_vector3f v[] ) const; 

int coordinates_at_vertex_cell ( const FEL_vertex_cell& c, 

FEL_vector3f * v) const; 

One can also convert between a structured mesh position s and physical coordinates 
via the call: 

int coordinates_at_structured_pos (const FEL_structured_pos& 

FEL_vector3f * v) const; 

Since the structured position includes a zone number, the “at structured pos” call can 
also be made on a multi-zone mesh, as long as the specified zone is structured. If the 
call is made on a mesh that does not have structured behavior, then the return value of 
the call will not be equal to 1. 
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11.6 Cell geometric properties 

The member functions which answer queries regarding the geometric properties of 
cells are relatively self-explanatory: 

bool cell_has_collapsed_edge (const FEL_cell&) const; 
int volume_of_cell ( const FEL_cell&, double*) const; 
int centroid_of_cell ( const FEL_cell&, 

FEL_vector3 f * ) const ? 

int longest_edge_length_of_cell (const FEL_cell&, 

float*) const; 

int closest_vertex_of_cell (cpnst FEL_phys_pos&, 

const FEL_cell&, 

FEL — vertex_cell*) const; 

The cel l_has^col lapsed.edge call tests whether two vertices on the same edge 
have coordinates that are exactly equal; in other words there is no epsilon used in the 
floating-point coordinate equality test to allow for nearly-equal values. 


11.7 Cell incidence relationships 

Given a cell c, an application may require cells incident to c. Two distinct cells c and 
d are incident if c is the face of d or vice versa. For example, for a given triangle c in 
a mesh, one may need the vertices of c, or perhaps the tetrahedra for which c is a face. 
The incidence relationships among cells can be visualized with a graph. Figure 11.2 
illustrates the incidence relationships for a small mesh in R 2 . The graph to the right 
contains a node for each cell in the mesh to the left. The nodes are organized into rows, 
each row containing cells of a particular dimension. The rows are ordered by ascending 
dimensionality: higher rows signify higher-dimension cells. A mesh containing 3-cells 
would have one extra row at the top. The example queries from above can be seen as 
starting at a particular node and following paths downward or upward. For example, to 
get the vertices of a triangle c, one could start at the node representing c and follow all 
the paths downwards. Likewise, in a mesh containing tetrahedra, one could start at a 
node representing a triangle c and follow the 0, 1, or 2 paths upward, depending on the 
number of tetrahedra that have c as a face. FEL meshes support queries based on cell 
incidence relationships via the calls up.cells and down.cells: 

int up_cells ( const FEL_cell& c, int d, int max, 

int* n, FEL_cell rc [] , int = -1) const; 

int down_cel Is (const FEL_cell& c, int d, int* n, 

FEL_cell rc[]) const; 

Both calls take a cell c as the first argument, and the dimension d of the cells one 
wants in return as the second argument. Both calls write cells in the array rc. The 
number of cells produced is written into n. The up.cells call also takes an argument 
max, which specifies the maximum number of cells that the user wants back. The 
terms “up” and “down” can be thought of either as going up or down in dimension. 
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Figure 1 1 .2: A small example mesh and its incidence graph. 


or as going upward or downward in the incidence graph. In topology, the up and 
down operations are known as star and closure , respectively. Note that the concepts 
of incidence relationships and up and down calls are not specific to a particular type 
of mesh; thus algorithms written in terms of up-cells and down.cells have the 
potential of working with many types of meshes. 


11.8 Adjacent cells 

A concept related to the incidence relationships between cells is adjacency, also known 
as a neighbor relationship. The mesh member function adjacent-cells supports 
queries requesting the cells neighboring to a given cell. The function signature for 
adj acent.cells is: 

int adjacent_cells (const FEL_cellS:, int*, FEL_cell []); 

where c is the cell for which to produce adjacent cells for, and nac and ac get the 
number of adjacent cells and the cells themselves, respectively. The most frequent us- 
age of the adjacent-cells call is with a 3-cell argument. For example, given a 
hexahedron c from a hexahedral mesh, adjacent-cells would return the hexahe- 
dra which share a quadrilateral face with c. Li k§wise ? given a tetrahedron c from a 
mesh (either a tetrahedral mesh or a hexahedral mesh with simplicial decomposition 
turned on), adjacent.cells will return the tetrahedra which share a triangle face 
with c. The adjacent -cells call is handy for algorithms that work breadth-first, 
starting from a seed cell. For example, one could construct an isosurface incrementally, 
processing cells outward from an initial 3-cell. 

The concept of adjacency has a more formal definition. First, returning to the graph 
in Figure 1 1.2, imagine that each vertex is the face of a special (-l)-cell, i.e. that there 
is an extra row beneath the vertex (O-cell) row with one node, and arcs from each 
vertex to the (-l)-cell node. Furthermore, if the cells in the top row are fc-cells, then 
imagine an extra row above the A:-cells with a single ( k + l)-ce11 that every A: -cell is 
the face of. Given this augmented incidence relationship graph, one can define the 
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adjacent cells of a cell c via the up and down operations described above. Let S u d be 
the set of ceils produced by going up one dimension and then down one dimension, 
starting with c. Let Sdu be the set of cells produced by going down one dimension and 
then up one dimension, starting with c. The cells adjacent to a cell c are the cells in 
(S u d fl Sdu ) - c. (There are more efficient ways to implement adjacent,cells, 
the preceeding definition is useful for its generality). 


11.9 Cardinality 

Since meshes are finite collections of cells, one basic query that the user might make 
is a count of the cells represented by a mesh. FEL provides this functionality via the 
method card: 

int card(int k) const; 

Given an integer argument k, card returns a count of the Ar-cells in a mesh. For 
instance, card (0) returns the number of vertices in a mesh. The value returned by 
card depends on whether simplicial decomposition is turned on or off. For example, if 
card { 2 } returns n when called on a structured hexahedral mesh with decomposition 
off, then card ( 2 ) will return 2n when decomposition is turned on. 


11.10 Cells and canonical enumeration 

For the Jt-cells in a mesh, one can imagine assigning a numbering so that each Ar-cell 
has a unique integer identifier. Such an enumeration could be handy, for example, 
for representing sets of A:-cells, since each integer identity number would consume 
less memory than an FEL.cell object. FEL supports a canonical enumeration of 
fc-cells, but not for every value of k and for every mesh type. The following chapters 
on structured and unstructured meshes list which enumerations are currently supported. 
Every enumeration follows the convention of going from 0 to card (k ) - i for Ar-cells. 
To convert between the integer representation and the FEL_cell representation, and 
vice versa, FEL provides the methods: 


int int_to_cell ( int i, int k, FEL_cell* c, 
int_ s = -I) const; 

int cell_to_int (const FEL ^ce l 1 & c , int * i ) const; 


The conversion from integer to cell is influenced by the simplicial decomposition mode 
currently set for the mesh. One can override the prevailing decomposition mode for the 
duration of the int_to.cell call by providing an optional final argument specifying 
a decomposition mode. For cell-to_int, the decomposition mode is inferred from 
the incoming cell type. 
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11.11 PLOT3D IBLANK 

The classes in the FEL mesh hierarchy, and the interface specified at the top hierarchy, 
are intended, as much as possible, to be independent of any particular mesh or file 
format standard. One case where FEL favors a particular standard is in its support of 
IBLANK. IBLANK is a standard defined by PLOT3D [WBPE92] where an integer 
“IBLANK” value is associated with each vertex in a mesh. Since meshes in FEL are 
currently all vertex-centered, having an integer at each vertex is equivalent to having 
an integer associated with each node in a field. The IBLANK value can serve one 
of three purposes. The first is as a flag indicating that the node associated with a 
vertex is invalid (indicated by an IBLANK of 0). The second IBLANK use is as a 
hint about overlapping zones in a multi-zone mesh. To indicate that a vertex may be 
within an overlapping zone z t PLOT3D specifies that the IBLANK value should be —z. 

(PLOT3D follows the FORTRAN style of numbering where the zones go from 1 to n 
rather than 0 to n - 1). The third usage of IBLANK values are to flag vertices that are 
on an impenetrable surface, signified by the value 2. An IBLANK of 1 is the default, 
signifying that the node is OK and that there is no overlapping zone information. FEL 
meshes provide access to IBLANK values via the calls: 

int iblank_at_cell (const FEL_cell& c, int i[]) const; 

int iblank_at_vertex_cell (const FEL_vertex_cell& c, 

int* i) const; 

int combined_iblank_at_cell (const FEL_cell& c, int* ci) const; 

int coordinates_ancL_iblank_at__cell (const FEL_cell& c, 

FEL_vector3f_and_int [ ] ci) const; 

int coordinates_and_iblank_at_vertex_cell {const FEL_vertex_cell&: c, 

FEL_vec t or3 f _and_int * ci) const; 

For a cell c, iblank_at_cell and iblank.at_vertex.cell produce IBLANK 
values, one for each vertex in c. The call combined_iblank-at_cell produces a 
single integer ci that is a bitwise combination of the following 4 flags: 

• FEL_PLOT3D_HASJBLANK_l 

• FEL_PLOT3D_HAS JBLANK.2 

• FEL_PLOT3D_HASJBLANK_0 

• FEL_PLOT3D_HASJBLANK_LT_0 

The combined IBLANK call is handy for quickly determining whether a more thorough 
analysis of the IBLANK values for a cell is necessary. The flags are designed so that 
a bitwise combination of them will result in an integer value between 1 and 15. FEL 
meshes also provide the “coordinates and iblank” calls, where one can request both 
types of data simultaneously. 

One can configure a mesh m so that a combined IBLANK value is returned by point 
location routines, using the call: 

m- > s e t ( FEL_RETURN_I BLANK , 1 ) ; 
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If locate finds a cell containing the given point, then the combined IBLANK for 
the cell is returned. Note that the combined IBLANK is an integer between 1 and 15. 
The locate function can still return other values, for example to indicate that the 
given point is outside the mesh. See FEL^re turns . h for a complete list of the return 
values. 


Chapter 12 

Structured Meshes 


In FEL, structured meshes consist of hexahedral cells (and all their faces), with a reg- 
ular organization such that the vertices of the mesh can be indexed by 3 indices, just 
as one would index a 3-dimensional array. The indices are usually written i, j, and fc. 
Topological information, such as incidence relationships among cells, is represented 
implicitly. Geometric information, such as the coordinates of vertices, can be rep- 
resented implicitly or explicitly, depending on the particular subclass of structured 
mesh. Subclasses of FEL.structured_mesh, in particular the curvilinear mesh 
subclasses, are some of the most heavily used classes in FEL. 


12.1 Simplicial decomposition 

When working with structured meshes, the user has the choice of 3 simplicial decom- 
position modes. Mode 0 corresponds to no decomposition. Modes 1 and 2 specify 
decompositions where each hexahedron is broken into 5 tetrahedra. There are two 5- 
tetrahedra decompositions possible for a h exah edron. In order for the decompositions 
to be consistent between each pair of adjacent hexahedra, the decomposition for each 
hexahedron must be the opposite of its adjacent neighbors. Thus, given the decompo- 
sition choice for one hexahedron in a structured mesh, the choices for all the remaining 
hexahedra are forced. FEL organizes the 2 decompositions in terms of “odd” and 
“even” vertices, where the odd and even designations come from the vertex i, j, and k 
indices. A vertex is even if the sum (i 4- j + k) is even, otherwise it is odd. In decom- 
position mode 1, the diagonals added to decompose the quadrilaterals go between even 
vertices. The decomposition choices for the 6 quadrilateral faces of a hexahedron leave 
only one possible tetrahedral decomposition. In decomposition mode 2, the diagonals 
go between odd vertices, and the hexahedral decomposition follows suit. 

12.2 Cell incidence relationships 

FEL-StructurecLmesh implements the structured mesh support for the up.cel Is 
and down.cells member functions. Table 12.1 summarizes the combinations of 
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Table 12.1: The supported combinations forup-cells and down_.ce 11s. Each box 
in the grid represents a combination of From cell dimension and To cell dimension. 
Combinations marked with d are supported by down.cells ( ) , combinations marked 
with u are supported by up_cells ( ) . The boxes marked • are trivially supported by 
either down-cel Is { ) or up.cells { ) . Empty boxes indicate unsupported combi- 
nations. 
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“from” and “to” cell dimensions supported. From the table, for example, one can see 
that since there is a d in row 3, column 0, one can use the method down.cells to get 
from a hexahedron (3-ceII) to its vertices (0-cells). The same support is available when 
simplicial decomposition is turned on; thus, given a subtetrahedron, one can get its 
vertices via down.cells. In general, the same “from” and “to” pairs are supported, 
regardless of the simplicial decomposition mode. 

The FEL structured mesh class also implements the adjacent cells method for struc- 
tured meshes. Currently adjacent-cells is implemented for hexahedra and subte- 
trahedra only. 


12.3 Canonical cell enumeration 

FEL.structured-inesh currently supports a canonical enumeration of vertices, 
hexahedra, and subtetrahedra. The enumeration for each type of cell corresponds to 
the ordering in which an iterator would produce the cells if iterating over the whole 
mesh. In terms of the cell z, j, and k indices, the i index would vary the fastest, k the 
slowest. If simplicial decomposition is turned on, then all the tetrahedra resulting from 
decomposing a particular hexahedron are numbered consecutively. 


12.4 Computational space support 71 ^ 

int coordinates_at_structured_pos (const FEL_structured_pos& s, 

FEL_vector3f* c) const; 

int jacobian_at_vertex_cell (const FEL_vertex_cell& v, 

FEL_matrix33f * m) const; 

int contravariant_phys_to_comp_vec tor (const FEL_vertex_cell& vc, 

const FEL_vector3f& pv, 
FEL_vector3f * cv) const; 

int contravariant_phys_to_comp_vec tor (const FEL_structured_pos& 

const FEL_vector3f& pv. 
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FEL_vector3f * cv) const; 


12.5 Structured mesh dimensions 

FEL provides access to the dimensions of a structured mesh via the member function 
get-Structured_dimensions ( ) : 

int get_structured_dimensions (int d[]) const; 

The function is defined as part of the interface of FEL_mesh, thus one does not have to 
cast an FELjnesh-ptr down to an FEL.structure_mesh_ptr in order to make 
the call. The function returns 1 on success or 0 if the call is inappropriate, e.g., if the 
call is made on an unstructured mesh. 

The get-Structured_dimensions ( ) call can also be useful for distinguish- 
ing objects that have structured mesh behavior from those that do not. We say “has 
structured mesh behavior” rather than “is a structured mesh” because not all meshes 
with structured mesh behavior are subclasses of FEL.structured_mesh. In par- 
ticular, a transformed mesh (Chapter 14) built in terms of a structured mesh has 
structured behavior, e.g., one can request the structured dimensions, yet it is not a 
structured mesh. Thus the call get_structured-dimensions ( ) is preferable to 
is_structured-inesh ( ) , since the fact a mesh is transformed should be transpar- 
ent to a routine that requires an object with structured mesh behavior. 


12.6 Axis-aligned structured meshes 

Axis-aligned structured meshes are meshes where the cells are aligned with the coor- 
dinate axes in physical space. Axis-aligned meshes in FEL may have either regular 
or irregular axes. A regular axis is an axis where the spacing between neighboring 
vertices is constant. An irregular axis is an axis where the spacing is not obliged to be 
the same throughout. An axis-aligned structured mesh where all the axes are regular 
is also known as a regular mesh . Axis-aligned meshes can be thought of as meshes 
defined by the Cartesian product of regular or irregular axes aligned with the physical 
space axes. 

Axis-aligned meshes have the advantage of requiring far less memory than curvi- 
linear meshes (described later in this chapter), since coordinates can be represented im- 
plicitly. (PLOT3D IBLANK values are assumed to be 1 for every vertex.) Axis-aligned 
meshes also have the advantage of more efficient point location, since the regularity of 
the geometry admits significant optimizations over the more general curvilinear mesh 
case. 

FEL provides classes representing axis-aligned meshes with all regular axes, 
meshes with regular x and y axes but an irregular 2 axis, and meshes where the z 
axis has dimension 1. The class constructors look like: 

FEL_regular_mesh ( int dO, int dl, int d2, 

char* nm = "regular^mesh" ) ; 
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FEL_regular_mesh(int dO, int dl, int d2, 

float sO, float si, float s2, 
char* run = R regular_mesh" } ; 
FEL_regular_mesh ( int dO, int dl, int d2, 

float sO, float si, float s2, 
float oO, float ol, float o2, 
char* run = "regular_mesh" ) ; 

FEL_regular_xy_irregular_z_mesh { int dO, int dl, int d2 , 

float oO, float sO, 
float ol, float si, 
float* coordinates2 , 

char* run = " regular_xy_irregular_z_mesh" ) ; 

FEL_regular_mesh2 (int dO , int dl, 

float sO, float si, 
char* = f, regular_mesh2 " ) ; 

The parameters dO, dl and d2 specify the dimensions of the mesh in I, J, and K. 

The run parameter gives the user the option of providing a character string name for the 
mesh. The second constructor for FEL_regular_mesh gives the user the opportunity 
to specify the spacing between adjacent vertices using the parameters sO, si, and s2. 

The origin for regular axes can also be specified using the arguments oO, ol, and o2 
in the third constructor. 

The class FELnregular_mesh2 is used for representing regularly gridded rect- 
angles. The arguments follow the same pattern as for the hexahedral meshes. 

An FEL«regularunesh2 instance lies in the z = 0 plane. The third com- 
ponent of FEL.phys.pos and FEL.structured.pos arguments is ignored by 
FEL-regular_mesh2. Thus, the locate member function, given a point p € R 3 , 
essentially projects p down to the z = 0 plane and then locates the quadrilateral to 
which p projects. 

The default interpolation mode with FEL_regularunesh2 is of type nearest 
neighbor. Simplicial decomposition is not currently supported for this class. 


12.7 Curvilinear meshes 

Curvilinear meshes are the most general type of structured mesh in FEL. The class 
FEL.curvi linear _mesh is actually an abstract class. The classes derived from 
curvilinear mesh that the user can instantiate are: 

• FEL.curvi linear_mesh_xyz -layout 

• FEL.curvi linear jnesh-xyzi.layout 

• FEL.curvi 1 inear _mesh.xy z i _f ielcLlayout 
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There are also paged curvilinear mesh classes, described in Chapter 23. The “lay- 
out” suffixes distinguish how the data describing the coordinates and in some 
cases IBLANK are provided. The xyz .layout mesh works with an array of type 
FEL.vec tor 3 f that contains the coordinates of the vertices. The x, y , and z compo- 
nents of each vertex are contiguous in memory, as they are in each FEL.vector3f . 
In terms of the structured mesh I, J, and K coordinates, the vertices are ordered so that I 
varies the fastest and K the slowest. IBLANK values are assumed to be 1 everywhere. 
The xyzi.layout mesh works with an array of type FEL.vector3f^and_int. 
The layout is the same as for the xyz mesh, except that each vertex has an 
IBLANK value interleaved with the coordinates. The xyzi.f ielcLlayout 
is constructed with a field where the node type is FEL.vec tor3f.and.int. 
When a FEL.curvilinearjnesh_xyzi_f ield.layout is queried for coordi- 
nates or IBLANK data, it in turn queries the field it was constructed with. The 
xyz i.f ield.layout curvilinear mesh is typically used to represent unsteady 
meshes; see Chapter 26. 

The constructors corresponding to the classes above are: 

FEL_curvilinear_mesh_xyz_layout ( 
int dO, int dl , int d2, 

FEL_vector3 f * xyz, 

const char* run = ,, curvilinear_mesh_xyz_layout " ) ; 

FEL_curvilinear_mesh_xyzi_layout { 
int dO, int dl , int d2 , 

FEL_vector3f_and_int* xyzi, 

const char* run = "curvilinear_mesh_xyzi^layout " ) ; 

FEL_curvilinear_mesh_xyzi_f ield_layout ( 
int dO, int dl, int d2 , 

FEL_vector3f_and_int_f ield_ptr xyzi_f ield, 

const char* run = n curvilinear_.mesh_xyzi_f ield_layout " ) 

Each constructor takes the I, J, and K dimensions of the mesh as the first 3 arguments. 
The next argument is specific to the particular data layout: the argument is either a 
pointer to a buffer or a field pointer. Thus if the user has a buffer in an appropriate 
layout, then he or she can construct a curvilinear mesh directly. 

In the future there may be other memory layouts supported by FEL via more sub- 
classes of FEL.curvilinearunesh. One layout that may be of particular interest 
is the case where the X, Y and Z coordinate values are in separate arrays, i.e., not in- 
terleaved together. Such a layout common in FORTRAN applications and in some file 
formats. 

The convention in FEL for mesh and field constructors is that any buffer provided 
as a constructor argument becomes the reponsibility of the FEL object to deallocate. 
Since FEL will use the destructor delete [ ] to do the deallocation, it is important 
that the user allocate the memory using the C++ allocation style for arrays, i.e., new 
[]. For instance, if the buffer is an array of FEL.vector3f, then the allocation 
should look something like: . 

FEL_vector3f * xyz = new FEL_vector3f [n_vertices] ; 
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If the user has a buffer allocated in some other manner, or in general if the user does not 
want FEL doing the buffer deallocation, then for a mesh m one can use the statement: 

m->set ( FEL_SUPPRESS_DEALLOCATION, 1) ; 

If the user specifies deallocation suppression, then he or she remains responsible for 
the management of the data buffers. See Chapter 7. 


12.8 Curvilinear mesh point location 

Curvilinear meshes inherit the interface of the two overloaded versions of locate 
declared in FEL_mesh. Both versions take an FEL_phys_pos as a first argument and 
a pointer to a cell where the result will be written as the last argument. The second 
version of locate takes an extra argument that is used as a start cell for walking 
point location. Given a start cell, locate for a curvilinear mesh will walk from cell 
to cell until a cell containing the given point is found. If locate with a start cell 
is unsuccessful, or if no start cell was provided, then curvilinear meshes uses four 
techniques to determine whether the mesh contains the given point. The techniques to 
locate a point p are: 

(1) Test if p is in the mesh bounding box. 

(2) Starting at the (computational coordinates) center of each of the 6 mesh boundary 
sides, use an adaptive vertex walk to get close to p. Starting from the walk 
destination closest to p, do a tetrahedral walk. 

(3) Starting from any of the destinations in (2) that have not been tried, do a tetrahe- 
dral walk. 

(4) For each 2-cell c on the boundary of the mesh, compute the centroid of c and the 
normal of c facing into the mesh. Choose the cell c from which p is visible and 
whose centroid is closest to p . Do a tetrahedral walk from there. 

The techniques are ordered by computational cost; the fourth technique in particular is 
relatively expensive. The user can control how much effort a curvilinear mesh m puts 
into point location using the call: 

m->set (FEL_LOCATE_EFFORT, level) ; 

where the mesh may use the level i technique if i < level. By default the level is 4, 
for multi-zone meshes the level for each zone is turned down to 3. 

The curvilinear mesh locate routines return 1 to signify success. If the mesh finds 
a cell containing the given point, it automatically checks the IBLANK values of the 
vertices of the containing cell. If any of the IBLANK values are 0, then the return 
value for the locate call is 0. 



12.9. CURVILINEAR SURFACE MESHES 


59 


12.9 Curvilinear surface meshes 

Another subclass of FEL_structured_mesh is FEL_curvilinear_sur- 
f acejnesh. The class represents surfaces in R 3 composed of quadrilaterals. The 
class has several constructors: 

FEL_curvilinear_surf ace_mesh ( int dO , int dl, 

FEL_vector3f * xyz, 
const char* = name_de fault) 
FEL_curvilinear_surf ace_mesh ( int dO, int dl, int d2 , 

FEL_vector3f * xyz, 
const char* = name_def ault) 
FEL_curvilinear_surface_mesh ( int dO, int dl, 

FEL_vector3f_and_int* xyzi, 
const char* = name_de fault) 
FEL_curvilinear_surf ace_mesh ( int dO, int dl, int d2 , 

FEL__vector3f_and_int* xyzi, 
const char* = name_de fault) 

The dO, dl, and d2 arguments specify the structured mesh dimensions. If three “d” 
arguments are provided, then exactly one of the three must be equal to 1, otherwise 
the d2 dimension is assumed to be 1. Following the “d” arguments is a pointer to 
a buffer with the coordinate data alone (FEL_vector3f *), or the coordinate data 
with IBLANK values (FEL_vector3f,and-int *). The user can query the mesh 
for coordinates and IBLANK values using the same “at” calls as for other types of 
meshes. The arguments should contain 0 in the “flat” dimension, i.e., the dimension 
given the value 1 in the constructor. 

Point location for surface meshes is defined to mean locating the 2-cell (i.e., quadri- 
lateral or triangle) whose centroid is closest to a given point location target p. If the 
distance from p to the closest cell centroid is shorter than the longest edge length of the 
cell, then the point location is considered successful. Currently it is not possible for the 
user to change the threshold for deciding whether a cell centroid is close enough. 

FEL also supports iterating over a surface. The one difference when working with 
a cell iterator and a surface is that the highest-dimension cells in the mesh are now 
quadrilaterals (or triangles if simplicial decomposition is turned on), rather than hex- 
ahedra. Otherwise, the iterators are initialized and used just as with other structured 
meshes. 
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Figure 12.1: A CFD visualization of a delta wing in a single-zone, structured mesh. 
The left side of the wing displays edges from the mesh, the right side show's contour 
lines tor a pressure derived field. See also Figures l.l and 13.1. (Data courtesy of Neal 
Chadorjiun, visualization courtesy of Tim Sandstrom. j 
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Unstructured Meshes 


An unstructured mesh can be thought of as a collection of cells, without the regular 
organization of a structured mesh. Currently FEL supports two types of unstructured 
meshes. The first, FEL.tetrahedral-mesh, contains tetrahedra, triangles, edges, 
and vertices. As part of the tetrahedral mesh construction, the user can provide sets of 
triangles defining particular surfaces. Constructing a “tetrahedral” mesh with sets of 
triangles but no tetrahedra is also allowed. The second type of unstructured mesh in 
FEL is an FEL.scattered..vertexjnesh. Scattered vertex meshes consist solely 
of vertex cells. 


13.1 Constructing a tetrahedral mesh 

The constructor for an FEL tetrahedral mesh looks like: 

FEL_tetrahedral_mesh ( int n_vert ices , 

FEL_vector3 f * coordinates , 
int n_special_triangles, 
FEL_vector3i* triangles , 
int* triangle_ids, 
int n_tetrahedra, 

FEL_vector4i* tetrahedra, 

const char* = " tetrahedral_mesh" ) ; 

Here n.vertices specifies the number of vertices in the mesh, coordi- 
nates is an array containing the physical space coordinates of each vertex, 
n.special-triangles specifies the number of items in the triangles array 
and in the triangle-ids array. Each FEL-vector3i in the triangles array 
specifies the three vertex indices of a triangle. The arguments n.tetrahedra and 
tetrahedra specify the number of tetrahedra, and the vertices of each tetrahedron, 
respectively. The vertex numbering in the triangles and tetrahedra arrays is 
expected to be FORTRAN style, i.e., the numbers should refer to vertices as if they 
were numbered 1 to n-vertices rather than 0 to n.vertices - 1. The buffers 
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Table 13.1: The combinations supported by FEL-tetrahedral .mesh for 

up-cells and down-cells calls. Each box in the grid represents a combina- 
tion of From cell dimension and To cell dimension. Combinations marked with d 
are supported by down.cells ( ) , combinations marked with u are supported by 
up.cells ( ) . The boxes marked • are trivially supported by either down.cells ( ) 
or up.cells ( ) . Empty boxes indicate unsupported combinations. 
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coordinates, triangles, triangle.ids, and tetrahedra passed in to the 
constructor become the responsibility of the tetrahedral mesh to deallocate when they 
are no longer needed. 


13.2 Cell incidence relationships 

Table 13.1 summarizes the currently implemented support for incidence relationship 
queries with tetrahedral meshes. From the table one can see, for example, that FEL 
tetrahedral meshes support queries requesting the vertices of a triangle (From 2 To 0), 
or the tetrahedra that a triangle is the face of (From 2 To 3). Queries for From/To 
combinations which are not currently supported return 0. 

The adjacent^cells call is supported for tetrahedral arguments only. 


13.3 Canonical cell enumeration 

FEL_tetrahedral_mesh currently supports a canonical enumeration of vertices, 
triangles, and tetrahedra. No enumeration is currently supported for edges. The enu- 
meration for each type of cell corresponds to the ordering in which an iterator would 
produce the cells when iterating over the whole mesh. 


13.4 Surfaces 

In structured hexahedral meshes, one can easily define a surface by holding one of 
the i, j , or k indices constant. Such surfaces are often used to represent entities such 
as the fuselage of an aircraft. Specifying a surface in a tetrahedral mesh is not as 
easy. To compensate for this, sets of triangles can be designated as representing a 
surface by the mesh generator, and a listing of the special triangles can then be included 
as part of the data set. For instance, the unstructured mesh file format defined by 
PLOT3D [WBPE92] supports the specification of triangle sets, each set with an integer 
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identity number. Note that the FEL_tetrahedral_mesh constructor starts with an 
array of triangles and an array of identity numbers, one for each triangle. As part of the 
construction, FEL-tetrahedral_mesh must determine the unique set of identity 
numbers and then rearrange the triangles into their corresponding sets. The user can 
query about the number of predefined triangle sets and get their identity numbers via 
the tetrahedral mesh calls: 

int get__n_triangle_sets { ) const; 

void get_triangle_set_ids { int ids[]) const; 

Using one of the predefined set identity numbers, one can initialize an iterator and 
loop over the triangles or vertices of a given surface. See the chapter on iterators 
(Chapter 16) for more details. 

13.5 Point location 

As with structured meshes, point location for a point p in tetrahedral mesh works by 
“walking” from 3-cell to 3-cell, until a cell containing p is found. The tetrahedral 
mesh locate is overloaded, just as in the structured mesh case, so that the user can 
provide a start cell for the walking search. If the walking search with a given start cell 
is unsuccessful, or if no start cell is given, then a global search over the whole mesh is 
done. 


13.6 Constructing a scattered vertex mesh 

The constructor for a scattered vertex mesh looks like: 

FEL_scattered_vertex_mesh ( int n_vertices , 

FEL_vector3f * coordinates, 
const char* = 

■ scattered_vertex_mesh" ) 

The array coordinates contains n.vertices of coordinates, i.e., coordinates for 
each vertex. The optional final argument gives the user the opportunity to give a specific 
name to the mesh. 

As with any mesh subclass, the scattered vertex mesh supports the standard mem- 
ber functions inherited from FEL_mesh. Cell incidence relationships are trivial, since 
there are only 0-cells. The canonical enumeration of the vertices is the same ordering 
as provided to the mesh constructor. Point location is defined as locating the ver- 
tex closest to the given search point. Interpolation is simply nearest neighbor. The 
FEL-cell-iter and FEL-vertex.ee 1 1-iter both have the same behavior, i.e., 
both iterate over vertex cells. 
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Figure 13.1: A CFD visualization of a fighter using an unstructured, tetrahedral mesh. 
The left side of the aircraft displays edges from the mesh, the right side shows contour 
lines for a pressure derived field. See also Figures 1.1 and 12. i. ( Data courtesy of 
NASA Lande\ Research Center, visualization courtesy of Tim Sandstrom ) 
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Transformed Meshes 


For domains where there is some type of symmetry, computational scientists often 
take advantage of the symmetry to model a smaller, fundamental domain. The results 
from modeling such a domain can then be replicated and transformed to fill the original 
domain. For example, in some turbomachinery studies, only one sector of a radially pe- 
riodic domain may be simulated. When visualizing the results from such simulations, 
the scientist may only need to visualize the results within the fundamental domain, or 
the scientist may wish to see the results replicated to look like the actual turbomachine. 
See, for example. Figure 14.1. In some cases, it may be sufficient to generate the 
graphics primitives for the fundamental domain and then draw the primitives repeat- 
edly, applying a different transformation each time. In other cases, replicating graphics 
primitives may not be enough. For instance, when using particle tracing visualization 
techniques, a scientist may be interested in seeing the traces continue beyond the fun- 
damental domaim In general, there are occasions when one would like to treat the 
simulation results as if they filled the whole original domain, without regard to any 
symmetry optimizations that may have been employed. 

FEL supports the representation of meshes with periodic symmetries through the 
subclasses of FEL.transf ormed_mesh. Transformed meshes are constructed with 
an original mesh and the data describing a particular transformation. Transformed 
meshes do not replicate the mesh data; thus one still enjoys most of the memory sav- 
ings due to not constructing a mesh for the whole original domain. At the same time, 
transformed meshes have the same interface as ordinary meshes and can be used just 
as non- transformed meshes, with no special treatment. The transformed mesh classes 
emulate rigid body transformations: translation and rotation. A transformed mesh can 
be constructed given any FEL mesh instance, including single-zone meshes, multi- 
zone meshes, and even other transformed meshes. See the mesh class hierarchy figure 
(Figure 1 1 . 1 ) for a refresher on the FEL mesh family. 
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Figure 14. i: A close-up of a periodic domain modeling a turbine [GBD961. The fun- 
damental domain used for the simulation would extend upward from the white regions 
at the base of the blades. Data courtesy of Karen Gundv-Burlet, visualization by Tim 
Sandstrom, 


14.1. HOW TRANSFORMED MESHES WORK 


67 


14.1 How transformed meshes work 

A subclass of FEL.transf ormedunesh is constructed with the mesh to be trans- 
formed and the data required for a particular transformation T. For mesh member 
functions that do not depend on the mesh geometry, such as card and up.cells, the 
transformed mesh simply delegates the call to the original mesh. For member functions 
involving the mesh geometry, such as calls producing coordinates or the bounding box, 
a transformed mesh makes the corresponding call on the original mesh and then applies 
T to the result. For point location, a mesh representing a transformation T applies the 
inverse transformation T -1 to the given point and then calls the locate routine on the 
original mesh, using the inverse transformed point as an argument. 

14.2 The transformed mesh subclasses 

The transformed mesh subclasses are illustrated under FEL-transformed-inesh in 
the mesh hierarchy (Figure 11.1). The constructors for transformed mesh classes are: 

FEL_translated_mesh(FEL_mesh_ptr m, const FEL_vector3f& t, 

char* run = " translated_mesh” ) ; 
FEL_rotated_mesh (FEL_mesh_ptr, const FEL_matrix33f & r, 
char* nm = " rotated_mesh” ) ; 
FEL_.x_rotated_mesh ( FEL_mesh_j?t:r m, float a, 

char* nm = H x_rotated_mesh"); 
FEL_y_rotated_mesh(FEL_mesh_ptr m, float a, 

char* nm = "y_rotated_mesh" ) ; 

FEL — z_rotated_mesh ( FEL_mesh_ptr m, float a, 

char* nm = " z_rotated_mesh H ) ; 

Each constructor takes a mesh m, and an optional name run. FEL.-translated.mesh 
emulates a translation by a vector t. The FEL_rotated_mesh classes emulate 
rotation transformations. The rotations are represented by a 3 x 3 matrix r such 
that given an original point p and a matrix R, the transformed point p ' = R * p* 
FEL also provides the classes FEL-x_rotated_mesh, FEL-y_rotated-mesh, and 
FEL.z.rotatedjnesh for representing rotations about the x, y, and z axes, respec- 
tively. The angle a should be in degrees. These classes make it easier to construct a 
mesh rotated about a particular axis, since one does not have to remember the matrix 
terms. The rotation subclasses also make it possible to do some optimizations, since 
the transformations can be computed in fewer floating-point operations than required 
for the general rotation case. 
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Chapter 15 

Multi-Zone Meshes 


Multi-zone meshes are represented in FEL by the class FEL_multi_mesh. A 
multi-zone mesh consists of 1 or more zones, where each zone is a subclass of 
FEL.single_mesh or perhaps a transformed version of some single mesh object 
(see Chapter 14). Note that the zones do not have to be all of the same type, for exam- 
ple, one could construct a multi-zone mesh containing both structured and unstructured 
zones. A multi-zone mesh with a zone that itself is a multi-zone mesh is not allowed, 
because objects such as cells that need to specify a zone contain a single integer for 
that purpose, and the integer is used to index into a single level of hierarchy. 

For most mesh member functions, a multi-zone mesh simply delegates the func- 
tion call to a particular zone. In particular, member functions that take an in- 
coming argument containing a zone number can be immediately delegated. The 
FEL classes FEL.cell and FEL-Structured-pos both contain zone integers, 
but FEL.phys.pos does not. By default, the zone in an FEL.cell or an 
FEL-S true tured.pos is set to FEL-ZONE-UNDEFINED. It is an error to call a 
multi-zone mesh member function with an incoming argument containing an unde- 
fined zone. 


15.1 Point location, IBLANK, and PLOT3D 

Point location is the most difficult task that a multi-zone mesh must support. Given a 
point p to locate, a multi-zone mesh must find a zone and a cell containing p. In gen- 
eral, the task is complicated by the fact that zones can overlap; thus p may be located in 
not just a single zone, but possibly multiple zones. When there is more than one zone 
containing p, then a “best” zone must be chosen. Thus, to be completely thorough, 
a multi-zone mesh would attempt to locate p in every zone and then choose the best 
result. The test-every-zone approach is effective, but too expensive for most applica- 
tions. FEL uses two strategies to improve upon the typical multi-zone point location 
performance, so that in most cases the test-every-zone approach is not necessary. 

The first strategy for accelerating point location is the same as that used with other 
FEL mesh types: provide a start cell for the walking point location. Since FEL cells 


69 


70 


CHAPTER 15. MULTI-ZONE MESHES 


contain a zone number, an initial cell argument effectively provides both a zone z to 
search first and a cell in z to start from. If the start cell is close to the location of the new 
point p, then it is possible that a cell containing p can be found quickly. Unfortunately, 
with a multi-zone mesh there is still the problem that the point p could be in more than 
one zone; thus even if the given start cell leads quickly to a cell c containing p, there is 
no assurance that c is the only cell containing p, or that c is the best cell containing p. 

The second strategy for accelerating point location addresses the overlapping zones 
issue, and in general the issue of which zone to try next if point location in a particular 
zone fails. The strategy relies on PLOT3D [WBPE92] IBLANK values. Recall that in 
PLOT3D, an integer IBLANK of -2 can be used to suggest another zone z to search. If 
point location in a paricular zone fails, then a negative IBLANK still makes it possible 
to avoid resorting to the test-every-zone strategy. Furthermore, if point location in a 
particular zone is successful, then the lack of negative IBLANK values can be taken to 
mean that there are no other zones overlapping a given point, thus no more searching 
is needed. 

IBLANK is also used by PLOT3D to flag nodes where the field data values are not 
valid (IBLANK of 0). As with curvilinear meshes, FEL returns an unsuccessful point 
location return code rather than a cell with a 0 IBLANK. In cases where there is more 
than one 0-IBLANK-free cell containing a given point, FEL must choose the best one. 
FEL chooses the cell with the smallest volume, since smaller cells typically come from 
higher-resolution meshes. • " ■■ 

FEL is designed to be explicit about features that are specific to a particular file 
format or CFD standard, as much so as is practical. For multi-zone mesh point lo- 
cation, FEL depends upon IBLANKs, thus the point location code is PLOT3D spe- 
cific. The PLOT3D dependency is made explicit in the FEL mesh class hierarchy: 
FEL_multi jnesh inherits from FELjnesh, and FEL_plot3djnulti jnesh in- 
herits from FEL_multi jnesh (see Figure 11.1). FELjnulti jnesh contains the 
vast majority of the interface and implementation for multi-zone meshes, the plot 3d 
class implements point location. In the future, other multi-zone mesh classes can be de- 
rived from FELjnulti jnesh, each with its own system for representing information 
analogous to IBLANK. 

15.2 Constructing a multi-zone mesh 

The class FELjnulti jnesh is an abstract class; to make a concrete multi-zone mesh 
object, one must currently use the PLOT3D class. The constructor looks like: 

FEL_plot3d_multi_mesh(int n_meshes, FEL_mesh_ptr* meshes, 

char* run = "FEL_plot3d_multi_mesh* ) 

The parameter runeshes specifies the number of meshes, meshes is an array of 
pointers to meshes. It is the responsibility of the multi-mesh class to deallocate the 
mesh pointer array when the multi-mesh is destructed. 

All FEL meshes return IBLANK values of 1 by default if no IBLANK data are 
provided. Thus one can construct an FEL_plot3djnulti jnesh instance even if the 
individual zones do not contain explicit IBLANK information. Point location efforts 
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will tend to be slower, since there will be no zone jumping hints, and situations where 
a point is contained in more than one zone will not be detected. 
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Chapter 16 

Iterators 


Iterators provide an interface for looping over the cells in a mesh. The iterator in- 
terface is independent of the particular mesh type, for example whether the mesh 
is structured or unstructured. FEL iterators work with the generalized cell object 
FEL.cell; thus iterators can represent not just hexahedra or tetrahedra, but also 
lower-dimensional cells such as vertices, triangles, and quadrilaterals. The class 
FEL.cell.iter provides the general cell iterator functionality. FEL also provides 
the class FEL.vertex.cell.iter, which is used in a manner very similar to 
FEL.cell.iter, except that the iterator only represents vertices (“vertex cells”). 


16.1 Basic iterator usage 

The basic use of iterators involves four operations: initialization, a done test, advancing 
to the next element, and dereferencing. The following code fragment illustrates all four 
operations: 

FEL_mesh_ptr mesh; 

FEL_vertex_cell_iter iter; 

FEL_vector3f c; 

int res; 

for (mesh->begin (&iter ) ; ! iter . done () ; ++iter) { 

res = mesh->coordinates_at_vertex_cell ( *iter , &c) ; 
cout << *iter << " coordinates: " << c << endl; 

} 

The initialization is handled by the begin method supported by all FEL mesh 
classes. The test whether the iterator is done is accomplished via the done ( ) method 
of iterators. Advancing the iterator is accomplished via the ++i ter call, and the itera- 
tor is dereferenced using the * iter syntax. Note how the dereferenced iterator can be 
used just as one would use an FEL.vertex.ceII argument, e.g., as an argument to coor- 
dinates.at.vertex.cell { ) . One can also call methods on *iter that belong 
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to the FEL.cell class, but one must be a bit careful about syntax. In C++ (and C), 
the operator has higher precedence than Thus, for example, the expression 
“*iter ,get_i jk ( ) M would parse the same as 44 * ( iter . get-i jk ( ) )”, which 
would not compile since the class FEL_cell_iter has no method get_ijk(). 
To access the methods belonging to the cell that an iterator represents, such as 
get.i jk ( ) , one should write 44 ( * i ter ) . ge t J. j k ( ) ’\ (The operator required to 
support -> syntax is not currently defined in the library; therefore one cannot write 
“iter->get_i jk ( )”) 

In a second example, we highlight two variations in the style of iterator usage: 


FEL_f ield_ptr field; 
FEL_vertex_.ce ll_iter iter, end; 
f ield->begin (&iter) ; 
field- >end j&end) -V' 1 

for ( ; iter != end; ++iter) { 

} " - : 1 * ' 


The first variation shows how one tests whether an iterator is done: rather than calling 
the done ( ) method, one can create a separate iterator object, initialized by the end ( ) 
method of meshes, and then compare the original iterator with the end object using 
the * = operator. The prefered FEL style is the former, i.e., use the done ( ) member 
function rather than making a separate end object. The latter style is provided for use 
in the future with the Standard Template Library (STL) [MS96]. The example also 
highlights the fact that the iterator initialization routines (begin ( ) and end ( ) ) are 
accessible via the field interface, so one does not have to get the mesh associated with 
a field in order to initialize iterators. 

The usage of FEL.cell-iter iterators is nearly identical to the vertex cell 
iterator above. Where FEL-vertex_cell-iter appears above, one would in- 
stead write FEL.cellJ.ter. Calls taking FEL_vertex_cell arguments would 
have to be replaced with the appropriate calls taking an FEL-cell argument. 
FEL-vertex_cell-iter objects always iterate over vertices. FEL.cell-iter 
iterators, on the other hand, by default, loop over the highest-dimension cell 
type in the mesh. The cell type produced by FEL.cel lJLter iterators is also 
a function of the simplicial decomposition mode for the mesh. So, for exam- 
ple, given a structured mesh containing hexahedra, an FEL.cel 1-iter produces 
hexahedra, or tetrahedra if simplicial decomposition is turned on. For a struc- 
tured surface mesh (FEL_curvilinear^urf acejnesh), FEL_cell-iter pro- 
duces quadrilaterals, or triangles if simplicial decomposition is on. Iterators over 
FEL.tetrahedral jnesh instances produce tetrahedra, regardless of the simplicial 
decomposition state. Finally, cell jteratorsover FEL-multijnesh instances produce 
the highest-dimension cell for each zone in die mesh. 
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16.2 Iterators and ordering 

FEL iterators allow the user to specify subsets of cells over which to iterate, but the 
order in which the user sees the cells is fixed. For structured meshes, iterators advance 
the I index the fastest, followed by J, then K. See Chapter 12 for a description of how the 
data members of an FEL.cell are set in order to specify a particular structured mesh 
cell. With unstructured meshes, only the I index is used. For the FEL-vertex-cell 
iterator on an unstructured-mesh m, the I index goes from 0 to m->card { 0 ) - 1. 
For the FEL.cell iterator on an unstructured_mesh m, the I index goes from 0 to 
m->card(3) - 1. 

In the case of multi-zone meshes, FEL iterators process the zones in ascending 
order, using the default order of the mesh for each zone to control how the cells are 
produced. 


16.3 Iterating over mesh subsets 

In some cases, one may desire to iterate over a subset of the cells of a mesh rather than 
every one. In order to describe a subset of a mesh, one typically must know about the 
mesh type, e.g., whether the mesh is structured or unstructured. Unfortunately, this 
implies that one must give up some of the the mesh type independence afforded by the 
default behavior of FEL iterators. Nevertheless, mesh subset iteration is important in 
certain cases. FEL supports mesh subset iteration via optional pairs of keyword- value 
arguments provided to the begin { ) initialization statement. 

Table 16.1 lists the keywords supported in iterator initialization. Also included in 
each of the S, U, and MM categories would be subclasses of the corresponding mesh 
class, for instance, FEL.curvi linear-mesh and FEL_regular_mesh would be 
included under the S category. The FEL_I_MAX parameter defaults to dim [ 0 ] - 1 
for a vertex iterator with a structured mesh. For unstructured meshes, FEL_I_MAX 
defaults to card ( 0 ) - 1. FEL_J_MAX and FEL JC.MAX are both set to 0. 

The following excerpt illustrates an example where the iterator produces the ver- 
tices on the K = 0 surface of a structured mesh, with a stride of 2 in the I and J dimen- 
sions: 

FEL_f ield_ptr field; 

FEL_vert_pos_iter iter; 

int res; 

res = f ield->begin ( &iter , FEL_K, 0/ 

FEL_I_STRIDE, 2, FEL_J_STRIDE , 2, 0); 

if (res 1=1) ... 

for ( ? ! iter .done () ; ++iter) { 

} 

Note that the begin ( ) statement above now has a return value, so that the library 
has the opportunity to indicate that one has provided initialization arguments that are 
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Keyword 

Default 

Mesh Types 

FEL.I-MIN 

0 

S, U, MM* 

FEL.IJMAX 

dimCO] - 1 

S, U, MM* 

FEL_I -STRIDE 

1 

S, U, MM 

FEL-J-MIN 

0 

S, MM* 

FEL.JJMAX 

dim[l] - 1 

S, MM* 

FEL.J-STRIDE 

1 

S, MM 

FELJK-MIN 

0 

S, MM* 

FEL-K.MAX 

dim [ 2 ] - 1 

S, MM* 

FEL-K.STRIDE 

1 

S, MM 

FEL.I 

(none) 

S, MM* 

FEL-J 

(none) 

S, MM* 

FEL.K 

(none) 

S, MM* 

FEL.UNSTRUCTURED-SURFACE 

(none) 

U, MM* 

FEL-ZONE 

0 

MM 


Table 16.1: The FEL iterator initialization keywords. The characters S, 

U, and MM in the column Mesh Types stand for FEL.structured-mesh, 
FEL.unstructured-mesh, and FEL-multi-mesh. classes, respectively. The as- 
terisk following MM indicates that the corresponding keyword is legal only if the ini- 
tialization is for a specific multi-mesh zone. 


not valid. For instance, the structured mesh initialization routine detects when the value 
for a parameter such as FEL-J-MAX is out of range, and returns a value not equal to 1 
to indicate this error. Note also that the final argument to the begin ( ) statement must 
always be 0 in order to indicate that are no more keyword/value pairs to follow. 

In the case where one is initializing an iterator for a multi-zone mesh, some initial- 
ization keywords are allowed only if one also specifies that the iteration should be over 
a particular zone, using the FEL -ZONE keyword. In Table 16.1, the rows where MM 
is followed by an asterisk designate keywords which are allowed only in conjunction 
with multi-zone meshes if one is selecting a specific zone. This restriction is due to the 
fact that most parameters, such as FEL-I-MAX or FEL_UNSTRUCTURED.SURFACE, 
typically only make sense when applied to a particular zone. 

The initialization of FEL.cell-iter instances is done using the same set of key- 
words as for FEL.vertex.cell iterators. The default values for the parameters 
are the same as for vertex iterators, except that the FEL-*_MAX values for structured 
meshes are initialized to dim [ 0 ] - 2,dim[l] - 2,anddim[2] - 2 for the I, 

J, and K indices, respectively. 

For unstructured meshes, the FEL-I-MIN and FEL-I-MAX parameters allow the 
user to control the first and last cell in the sequence of cells produced. One can also 
specify a stride via FEL-I -STRIDE. Parameters controling J and K are ignored by 
unstructured meshes. 
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16.4 Iterating over surfaces 

FEL iterators can also return vertices or cells from a surface. One can 
have a surface either because the underlying mesh is a surface mesh (e.g., 
FEL.curvi linear .surf ace_mesh) or because the iterator initialization key- 
words imply a surface. For structured meshes, one can specify a surface by holding ei- 
ther the I, J, or K index at a constant value. The iterator initialization keywords FEL-I, 
FEL_J, and FEL_K are used for this purpose. For example, in the code fragment above, 
the iterator initialization specifies the constant K = 0 surface. With structured meshes, 
iterating over a surface will produce quadrilaterals, or triangles if a simplicial decom- 
position mode is on. When iterating over a structured surface within a hexahedral 
mesh, the decomposition of the quadrilaterals into triangles matches the tetrahedral 
decomposition of the hexahedra in the original mesh. 

With unstructured meshes, sometimes the input file specifies sets of triangles 
which are taken to be part of a surface. Each set of surface triangles has an as- 
sociated indentification number. The user can access these triangle sets via the 
FEL JJNSTRUCTURED-SURFACE keyword, followed by a triangle set ID number. 

16.5 Iterators and time 

When working with time-varying data, one typically needs to initialize the time compo- 
nent of the cell represented by an FEL.cel 1-iter or an FEL_vertex_cell-iter. 
Both types of iterators support the methods set-time ( ) , set_physical-time ( ) , 
set-computational-t ime ( ) , and set.time.step ( ) . Using the set methods, 
one can initialize time for an iterator just as one would for an FEL cell. By default, the 
time associated with a cell is undefined. Note that FEL iterators do not currently iterate 
over time, i.e., advancing an iterator (using + +) does not change the time set by set 
calls described immediately above. On the other hand, the user is free to use the set 
calls within the iteration to manually change the time value. 
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Chapter 17 


Fields 


17.1 Fields in general 

The field class hierarchy is the centerpiece of FEL. Conceptually, a field is a contin- 
uous region of space over which some physical quantity is continuously defined. At 
any given spatial location in the field, the associated physical quantity may vary with 
time, or it may be time invariant. In finite-difference computational simulations, field 
values are calculated only at discrete points, but the spatial and temporal sampling 
are chosen to yield a reasonable approximation to a continuous physical system. The 
FEL field class hierarchy attempts to provide a set of objects and methods (data types 
and functions) which allow the user to treat finite-difference simulation data abstractly 
as continuous fields, with minimal regard for the actual underlying discrete spatial and 
temporal representation. Moreover, the FEL field class hierarchy provides mechanisms 
for combining different fields, and for converting one field type to another, according to 
various mathematical or physical identities, so that as field values are retrieved they are 
operated on by various functions, and can thus be returned directly in a form suitable 
for a particular task. 


17.2 Fields in context 

The FEL field class hierarchy is rooted in FEL_ref erence.counted.object via 
FEL-mutex_ref erence.counted-object, so that fields can be managed using 
reference counting, and in a thread-safe manner. Reference counting can reduce unnec- 
essary memory usage, a particularly desirable feature in the case of fields, which can 
be quite large. FEL_ref erence,counted_obj ect also harbors a character string, 
which can be used to name any of its descendants, including fields. 

Because fields are reference counted objects, they should be created only on the 
heap (via the new operator), using FEL’s smart pointers as handles. More on instanti- 
ating fields later. 
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17.2.1 Typeless fields 

Derived directly from the reference counting classes is FEL.f ield, the top of the field 
lineage proper. Unlike all of the other field classes, which inherit from it, FEL.f ield 
is not templated, and it contains not just a common interface, but common code shared 
by all fields, regardless of the type of field data. This node-type-independent code 
mostly provides an interface to the field’s mesh, and to two types of iterators which 
traverse the mesh. There are also some functions for translating between physical and 
computational time, and for retrieving various global physical quantities associated 
with the field. Factoring out such code common to all fields, from the more parochial 
code of particular field types, not only embodies a clean conceptual separation, but it 
also may reduce redundant code generation by compilers that take an ’’all or nothing 
approach to template instantiation. 

From a client’s point of view, FEL.f ields are important because they are the 
common ancestor of all fields — hence a pointer to FEL.f ield (typedef’d as 
FEL.f ield.ptr) can refer to a field of any type. The FEL.f ield pointer thus 
provides the means by which one can store handles to potentially diverse field types in 
a single homogeneous container (such as an array), or construct functions to operate on 
arbitrary fields. This feature is invaluable in writing general purpose code. Of course, 
FEL.f ields have such a generic interface that not much useful can be done with them 
without knowing some basic facts about their actual instantiations: methods to obtain 
these facts are provided by the FEL.f ield is.*.f ield ( ) interface. 

17.2.2 Typed fields 

Derived from the generic FEL.f ield is the templatized FEL.typed.f ield<T>. 
FEL.typed.f ield<T> is parameterized by <T>, which is a placeholder for the type 
of data provided at each node. This parameterization allows a single body of code to 
support instantiation of fields of any type, so long as the type supports a few basic 
arithmetic operations (see Chapter 21). FEL.typed.f ield<T> specifies the type- 
dependent interface common to all fields - primarily the at.* ( ) calls, which pro- 
vide lazy evaluation of field values at given locations in space and time. For many 
applications, the at.* { ) calls are the heart of the FEL user interface. In addition, 
FEL.typed.f ield provides a convenience function for producing an “eagerly eval- 
uated” field: this function precalculates field values at each node, and stores them into 
memory. By eliminating potentially redundant calculations, eager evaluation may be a 
logical choice if one plans on making heavy use of a highly derived field. 

There are at present eleven immediate descendants of FEL.typed.f ield (see 
Figure 17.1). The most heavily used types will be described in more detail in following 
chapters, but for now here is a brief overview. 


FEL.core.f ield<T> 

A field whose node values reside in memory. Core fields are typically produced by 
reading a PLOT3D solution file from disk, by executing get.eager.f ield( ) , or 
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JFEL core fieid<T>| 


jFEL_tlme_vafylng_field<f>| — jFELJInrH?_seHesJteid<f>l — | FEL_fixed_i n terval Jiroe_serl « J leJd < TV] 


fjFEL mesh as fiekJ<T> 


jFEL_cachedJi«kf<T>l 


JFEL hash cached fleid<T>] 


|FEl_fiftld >1 — IF El_typed_fiekjkT >~ >r~ -- jFE L_consta n t_f ield <T>] 


1FEL derived fiekJ<TQ> oj 


IFELi multi fiekj<T> 


l1FEL_differential_operatorJ?eki<TO t FROM>¥1 


nFEL_paged_field<T>] 


1FEL touch counted field<T>| 


Figure 17.1: A portion of the FEL typed field hierarchy. The subclasses of 

FEL.derived.f ield<T> and FEL.dif ferential-operator.f ield<TO, 
FROM> do not appear in this diagram. See Figure 20.1 for the differential operator 
subclasses. 


by explicitly associating a mesh with a node buffer in a core field constructor. See 
Chapter 18. 

FEL_paged-f ield<T> 

This is a field whose node values are paged in from disk on demand — useful or even 
imperative for extremely large data sets. See Chapter 23. 

FEL-derived.f ield<TO> 

Under the lazy evaluation paradigm, derived fields are essentially filters which trans- 
form field data as part of their retrieval. Derived fields are always built on top of other 
fields, and while the derivation chains can be arbitrarily long, there must always be 
at least one field at the bottom which can retrieve or manufacture a field value “au- 
tonomously” to get the whole thing started. Derived fields are templatized by the type 
of data they produce. See Chapter 19. 

FEL-dif f erent ial .opera tor.f ieldl<TO, FROM> 

FEL-dif ferential.operator.f ield2<TO, FROM> 

These are specialized derived fields which apply the nabla operator to scalar or vector 
fields, producing other scalar or vector fields, using either first- or second-order approx- 
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inflations for the required derivatives. The differential operator fields are templaUzed 
by their input and output types. See Chapter 20. 


FEL_constant_f ield<T> 

Constant fields return the same (constant) value from all locations. This can be useful 
for certain applications which want to combine fields algebraically — for example, one 
might shift the frame of reference of a velocity field by adding a constant velocity vec- 
tor at all points. Note that the same effect can be achieved with derived field mapping 
functions. 


FELjneshuas.f ield<T> 

This field type in effect creates a vector field whose entry at each node is just the po- 
sition vector of the node, as given by the mesh. Thus this field type constitutes an 
adaptor , allowing one to access a mesh using the field interface. Using this strategy, 
for instance, one can implement physical-space cutting surfaces on a mesh by extract- 
ing isosurfaces from the associated FEL_mesh_as_f ield. The mesh-as-field also al- 
lows one to convert computational coordinates to physical coordinates simply by using 
FEL_mesh_as_f ield : : at.structured.pos ( ) , but FEL provides a specialized 
and more efficient way of doing this: coordinates-at_structured_pos ( } , 
callable on any mesh. 


FEL_time_varying_f ield<T> 

This field type and its subclasses provide the additional interface necessary to support 
time-varying data. See Chapter 25. 

FEL_multi_f ield<T> 

This field type is basically a container class capable of managing multiple fields. It will 
be used primarily for supporting transformed fields. 


FEL_touch_counted_f ield<T> 

This field type simply sits on top of another field and records usage statistics from that 
field. Such statistics can be very informative during application development or tuning, 
as they may guide decisions about deployment of lazy evaluation, eager evaluation, or 
paged fields. 

FEL.cachecLf ield<T> 

This field is built on top of another field, and caches node query results from that field. 
Field queries on revisited nodes can be fetched from the cache, potentially saving time 
by eliminating redundant calculations on a highly derived field. 
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FEL Jiash.cached.f ield<T> 

Similar to FEL_cached_f ield, but the cached data is stored in a dynamically cre- 
ated hash table, instead of in preallocated memory as in FEL.cached-f ield. The 
hashing scheme may save a lot of unnecessary storage space in sparsely accessed fields, 
at the cost of slightly higher cache retrieval times compared to FEL.cached.f ield 
(but still potentially faster than pure lazy evaluation, on a highly derived field). 


17.3 Fields in detail 

17.3.1 Every field has a mesh 

A key data member of every field is its mesh, which is represented in the field by 
a FELjnesh-ptr, Every instantiated field must contain precisely one valid mesh 
pointer. It is not possible to create a field without supplying a mesh of the same cardi- 
nality as the data buffer, since each node value must have an associated spatial location 
before it is possible to carry out such basic operations as point location or interpolation. 

In contrast, a given mesh can be included by any number of fields, including none 
at all. A mesh by itself can support many purely geometric operations, whether or not 
any further data are attached to its vertices; and a single mesh can provide the spatial 
organization for many types of data, whether or not those data reside in memory or are 
constructed on the fly. A core field and any derived fields it supports will always share 
the same mesh. 

Thus in FEL there is a considerable asymmetry between meshes and (non- 
positional) node data. The FEL paradigm depends on a one-to-many but nevertheless 
tight binding between a mesh and its node data. The main business of creating a core 
field is establishing this binding; any derived field merely inherits its input fields 1 mesh. 

The FEL_mesh_ptr inside a field is protected, which means you can’t access it 
directly. The access function get jmesh ( ) , which you can call on any field, returns a 
pointer to the field’s mesh. You can use this pointer to construct other fields, or to call 
any of the publicly available mesh functions. For example: 

void f oo ( FEL_f ield_ptr field) 

{ 

FEL_mesh_ptr mesh = f ield->get_mesh ( ) ; 

FEL_vector3f lo, hi; 

mesh->get_bounding_box ( &lo , &hi ) ; 

FEL_vector3f_f ield_ptr coord_field = 
new FEL_mesh_as_vector3f_f ield (mesh) ; 

FEL_cell_iter ci; 

FEL_vector3f pvec [FEL_CELL_MAX__NODES] ; 

for {mesh->begin (&ci) ; fci .done() ; ++ci) 

{ 

coord_f ield->at_cell ( *ci , pvec) ; 
isosurface (pvec , ...); 

} 
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} 

Several of the commonly called methods on FEL-mesh have also been put in the 
FEL.f ield interface, so you can call them directly on a field without first retrieving 
the field’s mesh pointer. In these cases, the field merely forwards your function call to 
its mesh. Functions in this category include: 

int card (int); 
int get_n_zones { ) ; 

void set (FEL_set_keyword__enum, int) ; 

int coordinates_at_cell (const FEL_cell&, FEL_vector3f []) 
int coordinates_at_vertex_cell (const FEL_vertex_cell&, 

FEL_vector3f *) ; 

int convert_time (const FEL_time&, 

FEL_time_representation_enum, 

FEL_tIme* ) const; 


and the iterator functions: 

void begin (FEL_vertex_cel l_i ter*) ; 

int begin (FEL_vertex_cel l_i ter* , int, ...); 

void end (FEL_vertex_cell_iter* ) ; 

void begin (FEL_cell_iter*) ; 

int begin ( FEL__cell_iter* , int, - . . ) ; 

void end (FEL_cel loiter*) ; 

See Chapter 1 1 for details of these functions. 

17.3.2 Theat.M) calls 

The at.* ( ) calls allow one to retrieve field values at arbitrary locations and times 
within the computational domain. Locations and times can be specified in either com- 
putational or physical coordinates, with the previously noted exception that fields based 
on unstructured meshes don’t (can’t!) support computational spatial dimensions. Since 
the at.* ( ) calls are declared to return field values of a specific type, they are intro- 
duced to the field interface in the templatized FEL-typecLf ield<T>, the common 
ancestor of all instantiable fields. ' - -i _ 

In addition to being parameterized by the type of field value being retrieved, as 
just mentioned, the at_* ( ) calls are named according to the type of the field location 
being queried. This scheme was adopted, rather than overloading on location type, for 
reasons involving C++ function hiding. FEL’s design dictates that various specialized 
field types redefine certain functions that are initially defined as high as possible in 
the class hierarchy. The at_* ( ) calls are virtual, so that this redefinition gives rise to 
polymorphism. However, in C++ it is impossible to selectively override only a subset 
of a group of overloaded functions: if even a single member of the group is overridden, 
the rest are effectively hidden. This function hiding can only be overcome by redefining 
the entire group of overloaded functions, and for those that are unchanged, providing 
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explicit redirection up the class hierarchy. This solution is inefficient — it may take 
multiple (redirected) virtual function calls to reach the actual method — and it also 
results in needlessly cluttered class declarations and definitions. 

FEL largely sidesteps these issues by renaming each of the at.* ( ) calls according 
to the type of location being queried. This permits the different varieties of at.* ( ) to 
be independently overridden, at the cost of a slightly more cumbersome function name. 

Renaming functions according to their argument types is known as “name mangling , 
and is the strategy employed by C++ compilers to distinguish between overloaded func- 
tions. FEL is “partially mangling” the at.* ( 5 interface, but hasn’t abandoned over- 
loading altogether: the at.phys.pos ( ) calls are overloaded with respect to other 
argument types. 

The partial mangling and overloading in the at_* { ) interface results in 8 distinct 
calls: 

int at_vertex_cell (const FEL_vertex_cell&, T*) ; 
int at_cell (const FEL_cell&, T[]); 

int at_cell_interpolant (const FEL_cell_interpolant& , T[]); 
int at_structured_pos (const FEL_structured_pos&, T*) ; 
int at_phys_pos (const FEL_phys_pos&, T*) ; 

int at_phys_pos (const FEL_phys_pos& , FEL_cell_interpolant* , T*) 
int at_phys_pos (const FEL_phys_pos& , FEL_cell_interpolant&, 
FEL_cell_interpolant*, T*); 
int at_phys_pos (const FEL_phys_pos& , 

const FEL_cell_interpolant&, T*) ; 

The first three of these calls retrieve field values directly from nodes (or sets 
of nodes), and do not involve any spatial interpolation, at.vertex.cell () re- 
turns field values from an individual vertex in the mesh, and at_cell() returns 
an array of field values from the group of mesh vertices making up a cell. The 
cell can be any type supported by FEL (see Table 5.1), and since one cell type is 
FELXELL.VERTEX, at.vertex.cell ( ) is really just a specialized version of 
the more general at.cell ( ) . Because no spatial interpolation is necessary to fetch 
field values at cell vertices, at.cell-interpolant ( ) simply disregards the inter- 
polant part of the celLinterpolant, and in all other respects is essentially the same as 
at.cell ( ) . It is provided as a convenience to the user. 

The calls at.structured.pos ( ) and at.phys.pos ( ) can retrieve field 
values at arbitrary locations, as specified in computational (structured meshes 
only!) or physical coordinates, respectively. The requested location may not co- 
incide with a mesh vertex, so some type of spatial interpolation is necessary. In 
at.structured.pos ( ) the interpolation is always performed in computational 
space, using the computational coordinates supplied with the query. This is tantamount 
to isoparametric interpolation, discussed in Chapter 9. With at.phys.pos ( ) , the 
spatial interpolation is done in either computational space or physical space, depend- 
ing on the current interpolation mode (see Chapter 1 1 for details about setting the inter- 
polation mode). In either case, if repeated and relatively localized at.phys.pos ( ) 
queries will be made, much of the “set-up” work associated with the interpolation can 
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be profitably saved and potentially reused. Previous point location and interpolation in- 
formation can be cached in an FEL.cell-interpolant object, and the overloaded 
versions of at-phys_pos ( ) support the reuse of these hard- won data. For this rea- 
son, there is a hierarchy of at.phys.pos { ) calls, and the most appropriate call for a 
given field query depends on how much prior knowledge is on hand: 

• at-phys-pos {const FEL-phys-pos& , T*) does global point location to 
find the mesh cell enclosing the queried physical position, builds an interpolant 
based on that cell’s geometry, uses the interpolant to produce an interpolated 
field value which it stores into T*, then deletes the interpolant. This is the most 
expensive and most extravagant at_phys_pos () call, and is generally used 
only when an isolated field value is required. 

• at-phys_pos (const FEL_phys-pos&, FEL-cell-interpolant * , 
T*} does global point location to find the mesh cell enclosing the queried 
physical position, builds an interpolant based on that cell’s geometry, uses the 
interpolant to produce an interpolated field value which it stores into T*, then 
stores the interpolant into FEL.cell-interpolant*, whence it is returned 
to the client. This is just as expensive as the previous call, but at least now we 
have the possibility of reusing the information stored in the cell interpolant! 

• at_phys_pos (const FEL.phys_pos& , FEL_cell-interpolant& , 
FEL.cell-interpolant* T*) does “local” point location, i.e. starting 
in the cell of the supplied FEL_cell_interpolant&. This will converge 
much more quickly than a global search, if the queried physical point is 
nearby; and there is the added bonus that the interpolant part of the supplied 
FEL.cell-interpolant can be reused, if the queried physical point falls in 
the supplied cell. In any case, the current cell and interpolant are returned__vi_a 
FEL_cell_interpolant for the next go around. This is a safe and efficient 
method if successive at.phys.pos { ) calls are highly spatially correlated, as 
when integrating a stream line, for instance. Be warned, however, that the local 
search can be much slower than a global search, if the desired physical point is 
far from the supplied cell. 

• at_phys.pos (const FEL.phys.pos& , const 
FEL-ceii-interpolant &, T*) doesn’t do point location at all, but 
simply assumes the provided FEL.cell-interpolant is appropriate and 
goes ahead with the interpolation. This is certainly the fastest at.phys-pos ( ) 
call, but risky if you’re not sure the FEL.cell-interpolant matches the 
FEL-phys_pos. If it doesn’t, the at call may fail, but even worse, it may well 
succeed, and silently produce garbage that may or may not be easily recognized 
as such. Be sure you know what you’re doing if you use this streamlined 
at-phys.pos variant. 

All the at_* ( ) calls apply to time- varying fields, in which case temporal interpola- 
tion may be necessary ( may because temporal interpolation is bypassed if the requested 
time falls squarely on a timestep). The positional classes in the various at-* ( ) calls 
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all have an entry for time, which is simply ignored in time-invariant fields. See Chap- 
ter 25 for a more detailed discussion of time- varying fields. 

The at_* ( ) calls store the field values they retrieve into locations specified by a 
user-supplied pointer (T*). The return value of the function itself signifies the outcome 
of the retrieval: FEL.OK for success, and other values resulting from various mishaps. 
If the at ( ) call fails, the returned field value slot (T*) will probably be unchanged, 
and the returned FEL.cell.interpolant will probably be bogus. Since the T* is 
frequently reused, it may still point to a valid looking value, and be unwittingly used. 
Reusing a bogus FEL.cel 1-interpolant may result in gross inefficiencies, further 
failed at calls, or even a core dump (typical scenario: failure to create an interpolant re- 
sults in it being set to NULL, and then one blindly uses it in an at_phys_pos ( const 
FEL-phys _pos& , const FEL-cell..interpolant&, T* ) call ... SIGSEGV 
lies this way). For these reasons, good coding style, and general peace of mind, one 
should always check the return values of the at-M ) calls. 


17.3.3 Iterating over fields 

As a convenience, FEL-f ield supports the same iterator interface as FELjnesh. In 
fact, iterator methods invoked on a field are merely forwarded to the field’s mesh. This 
shorthand somewhat compromises the abstraction of a field as a continuous domain, 
but allows a more relaxed coding style in which one simply makes all function calls off 
a field, without having to remember which ones are more mesh-related. For instance: 


int foo (FEL_f loat_f ield_ptr field) 

{ 

float f [FEL_CELL_MAX_NODES] ; 

f ield->set ( FEL_SIMPLICIAL_DECOMPOSITION, 0) ; // 
FEL_cell__iter ci; 

for ( f ield->begin (&ci ) ; ici.done; ++ci) // 

{ 

field- >at_cell ( *ci , f ) ; // 

if (any (f ) <0 .0) 

f ield->coordinates_at_cell ( *ci , c ) ; // 

f ield->at__phys_pos {average (c) ,8c f) // 

cout << (*ci) 

<< "has negative vertex: * 

<< "scalar at centroid = H 
« f 
« endl ; 

} 

} 


forwarded 

forwarded 

field method 

forwarded 
field method 


See Chapter 16 for a more complete discussion on the capabilities and uses of 
iterators. 


88 


CHAPTER 17. FIELDS 


173.4 Eager fields 

FEL defaults to “lazy evaluation” of field values. This means that all field values FEL 
returns via the at_* ( ) calls, and any intermediate values required for their evaluation, 
are generated only when needed to satisy a given at.* ( ) call. The sole exceptions to 
this rule are the FEL_core_f ield nodal values, which reside in a buffer in memory, 
or perhaps on disk in the case of a FEL_paged.f ield. The creation of a derived field 
merely sets up an appropriate filtering mechanism, which is not pressed into action 
until derived values are requested. Particularly in the case of highly derived fields, 
lazy evaluation saves a lot of storage space and setup time, at the expense of slower 
turnaround when any values are actually needed, and potential redundant calculations 
if the same locations are queried over and over again. 

Sometimes it might make better sense to evaluate derived fields at all their mesh 
vertices and store these results in memory. Then when node values are needed, FEL 
merely has to fetch them, for immediate return, or for derivative or interpolation pur- 
poses. This “eager evaluation” strategy essentially converts a derived field into a core 
field. It requires more storage space, but particularly in the case of highly derived 
fields, can provide faster immediate access, and can also provide longer term savings 
by eliminating redundant calculations if the same locations are repeatedly queried. An 
eager field is a field level cache. 

One creates an eager field by invoking get.eager.f ield { ) , which is 
a method on FEL_typed_f ield, so it can be called on any typed field, 
get .eager _f ield ( ) returns a FEL.po inter to the same type of field that calls it. 
Thus, assuming lazy.vector.f ield is already defined: 

FEL_vector3 f_f ield_jptr eager_vector_f ield = 
lazy_vector_f ield->get_eager_f ield { ) ; 

In the case of time-varying fields, get.eager.f ield ( ) requires an argument 
specifying the computational time of the newly produced field; see Chapter 25 for 
details. If get^eager.f ield ( ) fails — most likely because it can’t allocate enough 
memory or, in the time-varying case, if the time is invalid — it returns NULL. Be 
sure to check the return value of get_eager_f ield{ ) before you do something 
embarassing like trying to dereference a NULL pointer. 

Eagerly evaluated fields are probably a logical choice if you are faced with some 
combination of the following: ' 

1 . you need the quickest possible access to a field 

2. you are repeatedly accessing lots of locations in the field 

3. you have a very highly and expensively derived field 

4. you have enough memory 

An application supporting interactive sweeping of a gridplane through a velocity x 
vorticity magnitude field is a good example meeting the first three criteria above. 
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17.3.5 Field type “informant” functions 

The FEL field class hierarchy strives for polymorphic behavior, but sometimes you just 
can’t pretend anymore and you simply must know what specific type of field is really 
being represented by some generalized FEL field pointer. For such desperate times, 
FEL provides type-specific “informant” functions which return true or false depending 
on whether the invoking field meets the criterion encoded in the function name. There 
are six such functions: 

bool is_core_f ield( ) ; 
bool is_ time_series_f ield{ ) ; 
bool is_ f loat_field( ) ; 
bool is_vector3f_f ield( ) ; 
bool is_plot3d_q_f ield( ) ; 

bool varies_with_time ( ) ; 

These functions are methods on FEL.f ield, so they can be called on any type of field 
whatsoever. Note that a given field can return true for more than one of the queries; 
only the return values from is.f loat.f ield ( ) , is-vector3f_f ield ( ) , and 
is-.plot3d.q-f ield ( ) are mutually exclusive. 

An affirmative response to an is.*-f ield ( ) call indicates that the queried field 
really is an object of the queried type, so that a valid downcast to the queried type is 
possible. In fact, type confirmation before downcasting is the primary intended use of 
the is.*.f ield ( ) calls. 

On the other hand, varies.with.time ( ) can return true for any of the field 
types which can be built on top of a time-varying field. The varies.with.time ( ) 
query is generally used as part of a switch or optimization in a general application, since 
the methods and resources for managing steady and unsteady data are so different. 

Note that on a given field, varies-.with.time { ) can return true and 
is-time-series.f ield ( ) may return false. This would happen, for example, if 
the field in question were a derived field built on top (directly or indirectly) of a time se- 
ries field. Keep in mind that an is_*.f ield( ) call tells you where a field is declared 
in the FEL class hierarchy, whereas varies.wi th_t ime ( ) tells you only that a field 
has a time-varying member somewhere in its individual lineage (the group of fields that 
are “chained” together by successive definitions). In biological terms, is.* .field { ) 
calls are queries about phylogeny, and varies.wi th.t ime ( ) pertains to ontogeny . 

Using the informant functions, one can be tempted to write code like this: 


int f oo ( FEL_f ield_ptr field) 

{ 

if ( f ield->is_f loat__f ield ( ) ) 

do_scalar_stuf f (FEL_FIELD_TO_FLOAT_FIELD_CAST (field) ) ; 
else if (f ield->is_vector3f_f ield( ) ) 

do_vector_stuf f ( FEL_FIELD_TO_VECTOR3F_FIEIiD_CAST (field) ) 
else 

return 0; 
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return 1; 

} 


Suit yourself, but be aware that this explicit style may become hard to maintain if 
used with abandon. If you have lots of type-dependent code, it may be better to localize 
the branch points, say by pushing them all inside a “Factory” class, which manufactures 
polymorphic objects that can be passed to type-independent functions, 

17.3,6 Odds and ends 

Min/max values 

Sometimes it is useful to know the minimum and maximum values of a scalar field, for 
example if one is fine-tuning a transfer function. To this end, FEL provides the method 
get_min_max ( ) , which returns the minimum and maximum nodal values of a scalar 
field: 

int res ; 
float min, max; 

res = some_f loat_f ield->get_min_max (&min, &max) ? 
if (res == FEL_0K) ... 

For time- varying fields the method requires an additional argument specifying the 
computational time, and as usual, the function return value signals the final outcome of 
the call. Although the method is declared on FEL_typed_f ield, it produces mean- 
ingful results only on float fields, i.e. fields on which is-f loat.f ield( ) returns 
true. For all other field types, get jninunax ( ) prints a message on standard error, 
leaves the pointer arguments untouched, and returns something_other than FEL.OK. 

In addition to taking care of the dirty work of iterating over the field in search of ex- 
trema, get_min_max ( ) caches the minimum and maximum field values once it finds 
them. On any given field, the second and subsequent invocations of get juinunax 
merely conjure up the cached values. This is very fast, and makes it unnecessary for 
the client to store the minimum and maximum field values explicitly. 


User data 

All FEL fields have two slots for user-defined data: 

int user_type; 
void* client_data; 

These are part of the top level FEL_f ield interface, so they are available in any 
field. 

The user_type entry is a single user-managed integer that is meant to be used as 
a simple type-tag if, say, the user wants to categorize fields differently than FEL does. 
One scenario has the user assigning to the tag at field creation time, using a custom 
enumeration. Functions can then be written to branch in their treatment of the fields 
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after examining the tag (but see warnings above in Section 17.3.5. The tag could also 
be used to track individual instances of fields, rather than types . 

The client-data is a general purpose pointer that can be used to associate ar- 
bitrary data with a field. One simple possibility for using client-data is just a 
generalization of the user- type scenarios. The user could attach a struct to a 
field, which contained type and instance records, and perhaps some other useful de- 
scriptive information. Another idea is to use the client -data hook to turn tradi- 
tional data-driven frameworks inside-out: Instead of constructing a framework which 
manages fields , by explicitly associating them with particular meta-data, visualization 
techniques, and so on, one could stuff all the bookkeeping inside the field itself, to- 
gether with appropriate behaviors for interacting with certain environments — so that, 
in the extreme, fields could largely manage themselves. Developers should keep in 
mind, however, that while the cl ient.dat a mechanism allows fields to be arbi- 
trarily enhanced and extended, such modifications are likely to be relatively domain 
specific, and should not be confused with enhancements and extensions to the Field 
Encapsulation Library itself. 
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Chapter 18 

Core Fields 


The previous chapter discussed the generic interface declared for all fields on the 
FEL.f ield and FEL.typed.f ield base classes. This chapter presents a more spe- 
cific treatment of a particularly important subclass of typed field called a core field. 


18.1 What are core fields? 

Core fields (FEL-core.f ield<T>) are fields whose node values reside in memory. 
This is in contrast to the various sorts of derived fields, for which node values are 
generated only on demand. Core fields live at the roots of derivation chains (see 
Chapter 19) and provide raw values which derived fields modify. 

Core fields are produced by reading a solution file from disk, by executing 
FEL.typecLf ield: :getJeager_f ield{ ) , or by explicitly associating a mesh 
with a node buffer in a core field constructor. 

Reading solution files from disk is_ the subject of Chapter 22, and 
get-eager-field ( ) is discussed in Chapter 17. Here we will describe how to 
construct a core field “manually” in memory. This operation is essential, for example, 
if one wants to load data directly from a field simulation into FEL, or if one wants to 
write a file reader for FEL. 


18.2 The node buffer 

As with any other type of field, construction of a core field requires a preexisting mesh. 
Meshes may be created by reading in a file from disk (see Chapter 22), or by declaring 
a FEL .regular .mesh. Details for creating a mesh “manually” in memory are given 
in Chapter 12 (for structured meshes) and Chapter 13 (for unstructured meshes). 

Assuming a mesh is already on hand, the central task of creating a core field is 
allocating a buffer and filling it with the node values. In most cases, the buffer should 
be allocated on the free store using either the new operator or operator new (or 
operator new [ ] , if your compiler supports it). This is because in most cases, once 
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the core field is created, memory management of the node buffer is turned over to FEL, 
and when the reference count to a core field reaches zero, FEL will try to deallocate its 
node buffer using delete [ ] . Attempting to delete [ ] memory which hasn’t been 
obtained via new is disastrous. 

If you want to retain responsibility for memory management of core field node 
buffers, you can do so by requesting that the core field “suppress deallocation”. If this 
option is chosen, the node buffer is left untouched when the core field is destructed. 
In this case, therefore, the memory need not be allocated with new — the memory 
may be obtained by some other dynamic allocator, or it may be static. In any event, 
the user is responsible for maintaining a handle to the node buffer, and for freeing the 
memory if so desired. Of course, the user should not deallocate the node buffer while 
the core field is still in use! Directions for choosing the “suppress deallocation” option 
are given below. 

Core fields are templatized, and can be constructed with any node type T supporting 
a few basic arithmetic operations (see Chapter 21). Since the node buffer must hold 
a type T for every vertex of its associated mesh, the buffer must be at least of size (in 
bytes): 

mesh->card ( 0) * sizeof(T) 

The elements of the node buffer must be arranged in the same order as the corre- 
sponding vertices of the associated mesh. For structured meshes of dimension (idim, 
jdim, kdim), this means that the sequential memory layout has idim varying most 
rapidly, and kdim most slowly. C/C++ programmers should be wary here, as this 
column-major convention may seem out of place (but in general is more convenient 
for data coming from predominantly FORTRAN solvers). For unstructured meshes, 
the sequential layout is arbitrary and completely determined by the vertex ordering in 
the mesh. In both the structured and unstructured cases, FEL vertex cell iterators fol- 
low the linear arrangement of mesh vertices in memory, so if you are creating a node 
buffer whose values depend on vertex coordinates, an FEL.vertex.cell.iter can 
be employed to guarantee vertex -node correspondence. 

For core fields with multi-element node types, based on either structured or un- 
structured meshes, the multiple elements associated with a given node should be layed 
out as a contiguous group in memory. Note that this natural layout differs from certain 
file formats, including the PLOT3D solution file format, in which vector components 
are grouped together. 

18.3 Constructors and suppressed deallocation 

There are three core field constructors. With a preexisting mesh and an allocated and 
filled node buffer, one creates a core field with one of the following two constructors: 
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FEL_core_f ield (FEL_mesh_ptr , T*, char* = "core_f ield" ) ; 
FEL_core_f ield (FEL_mesh_ptr , T* , bool, char* = "core_f ield" ) 

The FEL-mesh_ptr argument takes a pointer to the intended mesh of the core field. 

T* should point to the head of a node buffer of type T which you have allocated and 
filled with type T's, according to the scheme outlined above. The bool argument in 
the second constructor is used to indicate whether or not you want to suppress deallo- 
cation of the node buffer when the core field’s reference count reaches zero: true for 
suppressed deallocation (you are responsible for deallocation), false for automatic 
deallocation (FEL assumes management of the node buffer). If you use the first con- 
structor, “suppressed deallocation” defaults to false. The last core field constructor 
argument is an optional name, which defaults to “coreJield” if you don’t supply a more 
imaginative name of your own. 

The third core field constructor — 

FEL_core_f ield(FEL_mesh_ptr, FEL_pointer< FEL_core_f ield<T> 
char* = " shared_core_f ield" ) ; 

— takes a pointer to a preexisting core field rather than a pointer to a node buffer, 
and creates a new core field by binding the node buffer of the preexisting core field 
to the incoming mesh represented by FELunesh.ptr. This results in two core fields 
which share a single node buffer but associate the data with different meshes. For this 
to make sense, the different meshes must have the same number of vertices (meshl- 
>card ( 0 ) == mesh2 ->card ( 0 ) ) — if this is not the case, bad things will hap- 
pen. In a “shared core field” the node buffer is a shared resource, so the suppress 
deallocation option defaults to true. 

As we have just seen, the suppress deallocation option on a core field is set at 
construction time, either implicitly or explicitly, but it can also be changed at any time 
by a set call, if you change your mind about node buffer memory management: 

// FEL deletes buffers: 

my_core_field->set (FEL_SUPPRESS_DEALLOCATION, 0) ; 

/ / you 7 re on your own : 

my_core_f ield->set ( FEL_SUPPRESS_DEALLOCATION, 1) ? 

It’s probably not a good idea to set suppress deallocation to false (0) on a shared core 
field. 


18.4 An example 

Here is a rather contrived example which constructs a core field “manually”. See 
Chapter 12 for a description of FELuregularunesh, and Chapter 16 for details on 
FEL_vertex..cell-iter. 
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#include "FEL.h" 

int main() 

{ 

int idim=3, jdim=4 , kdim=5 ; 

FEL_mesh_ptr mesh = 

new FEL_regular_mesh(idim, jdim, kdim, 1,1,1) ; 

FEL_vector3f * node_buffer = 

new FEL_vector3 f [mesh->card ( 0 ) ] ; 

assert (node_buf fer) ; 

int n=0; 

FEL_vertex_cell_iter iter; 

// fill node buffer: 

for (mesh->begin ( &i ter ) ; liter. done () ;++iter) 

{ 

int i = ( *iter) [0] ; 

int j = { *iter) [1] ? 

int k - ( *iter) [2] ; 

node_buf f er [n++ ] . set (atan2 ( j , i ) , 

atan2 (k, sqrt ( i*i + j * j ) ) , 
sqrt ( i*i+ j * j+k*k) ) ; 

} 


// create core field 

// (suppress deallocation default is false) 
FEL_vector3f_f ield_ptr my_core_f ield = 
new FEL_core_f ield<FEL_vector3f > 

(mesh, node_buffer) ; 

FEL_vector3f v; 

for (mesh->begin (Sciter) ; liter. done () ;++iter) 

{ 

my_core_f ield->at__vertex_cell ( *iter, &v) ; 
cout << *iter << " : " << v << endl; 

} 

my_core_f ield = NULL; // node_buffer delete [] 'd 
return 0; 
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18.5 get_eager_field() 

Core fields are also generated by the get .eager _f ield { ) call. This method, avail- 
able on any field, evaluates the field at every vertex and writes the results into a newly 
allocated node buffer, which is bound to the mesh of the calling field. This function is 
typically used to convert a lazily evaluated derived field into a core field, for quicker 
access time, but could conceivably be used to duplicate an existing core field, in the 
unlikely event this was useful. More verbiage on get_eager_f ield can be found in 
Chapter 17. 
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Chapter 19 


Derived Fields 


19.1 What is a derived field? 

A “derived field” is a field whose values are computed (derived) from one or more other 
fields. The derivation referred to here involves some kind of mathematical mapping, 
and should not be confused with the C++ sense of derivation (although the FEL derived 
field classes do inherit from other classes). 

A derived field is essentially a filter, built on top of some number of component 
fields . The derived field and all its component fields necessarily share the same mesh. 
When the derived field is queried for a value at some location, it forwards the request 
to its component fields, gathers the values returned by the component fields, and then 
transforms these field values according to some function before passing them back to 
the caller. 

The function that the derived field uses to produce its values from its component 
field values is called a mapping function. As its name suggests, the mapping function 
maps component field values into derived field values. The mapping functions are 
either predefined by the library (more on these below) or user-defined. User-defined 
mapping functions allow arbitrary transformations of field data to be encapsulated in 
the retrieval process, so that field queries on derived fields can return data directly in a 
form suitable for a particular task. 

The component fields of a derived field can be of any type, including other de- 
rived fields. This last possibility allows “derivation chains” of arbitrary depth, and at 
any level the three types of component fields can be freely mixed. Naturally, several 
different derived fields can share component fields. Thus an application can build an 
arbitarily complex set of derived fields, but their relationships can always be described 
by a directed acyclic graph of one or more components. Since queries on derived fields 
are always forwarded to component fields, any given derivation “lineage” must even- 
tually terminate in a field type capable of producing a value autonomously. The field 
types that fit this bill are core fields (which pull a value from memory), paged fields 
(which pull a value from memory or disk), “mesh as fields” (which use mesh coordi- 
nates as field values), and constant fields (which simply return the same value for all 
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queries). 

19.2 Lazy vs. eager evaluation 

By default, evaluation of derived field quantities is completely demand driven. That is, 
no field values are generated when a derived field is created: derived field quantities 
are computed only in response to a client query. This strategy is variously described as 
lazy evaluation , or deferred evaluation , or a pull model . 

Because of the lazy evaluation paradigm, creation of derived fields requires very 
little time and storage. Internal creation of a derived field consists mainly of registering 
the component fields and mapping function, which are represented by pointers in the 
derived field. On our workstations (SGI RIOKs), creation of a derived field requires 
about 30 microseconds, and about 200 bytes of memory. Thus it is quite reasonable 
for an application to present many derived fields (hundreds, or more!) to a user as 
selectable options, creating them either at startup time or on demand, even if most of 
the fields are never used. 

The downside of the lazy evaluation scheme is lengthier derived field value retrieval 
times, compared to core fields, since the derived field must manufacture its values on 
demand. For simple derived fields (single core component field, mapping function re- 
quiring only local values) retrieval times are about 25% longer than core field retrieval. 
As the derivation chain comprises more component fields and more involved mapping 
functions, the derived field retrieval times increase accordingly. 

If retrieval time is at a premium and ample storage space is available, it may be 
desirable to precalculate and store derived field values across the entire domain, or at 
least retain and reuse those derived field values which are generated to satisfy client re- 
quests. The first of these aims can be accomplished by way of get.eager.f ield ( ) , 
a method on FEL_typed_f ield which iterates over all vertices of its invoking field 
and writes their values into newly allocated storage. This converts a derived field 
into a core field 1 , and may be the logical choice if one wants to minimize retrieval 
time on a highly derived field. The second option above can be realized by creating a 
FEL.cached.f ield on top of a derived field, which retains derived field values in 
a cache, and attempts to satisfy client requests from the cache before deriving anew. 
Cached fields may be useful if a derived field is going to be accessed repeatedly. Ex- 
amples of creating eager and cached fields are given below. 


19.3 Mapping and interpolating 

There are two distinct types of derived fields in FEL: FELjnap-then.inter- 
polate-derived-field* and FEL-interpolate.then_map_der- 
ived.field*. As the names suggest, the difference between these two types 
depends on whether the mapping function is applied before or after any necessary 
interpolation. More specifically: To satisfy a general physical or computational 
space query (i.e., at.phys-pos { } or at_structured_pos ( ) ), a field must first 

* ... or can also be used to generate a copy of a core field, in the unlikely event this was desired. 
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perform point location to determine which gridcell contains the query point. Then 
field values are obtained for each vertex of this enclosing cell and used to interpolate a 
field value at the enclosed point. 

• In an FEL_map.then-interpolate-derived-f ield*, the derived field 
queries its component field(s) at each vertex of the enclosing cell, applies the 
mapping function at each vertex to convert component field quantities to derived 
field quantities, then interpolates the derived field quantity at the query location. 

• In an FEL.interpolate.then_map.derived_f ield*, the derived field 
queries its component field(s) at each vertex of the enclosing cell, interpolates 
the component field quantities at the query location, then applies the mapping 
function to convert the component field quantities to the derived field quantity. 

The relative order of mapping and interpolation can make a big difference 
in the final derived field value, particularly if the mapping function is nonlinear. 
Which ordering is most "appropriate” depends on the problem at hand. In ad- 
dition to numerical implications, there are efficiency considerations: mapping be- 
fore interpolation invokes the mapping function on every vertex (potentially expen- 
sive for a complex mapping function) but interpolates only a single variable; in- 
terpolating before mapping invokes interpolation on all component field values (po- 
tentially expensive if there are several component fields, with complex node types) 
but applies the mapping function only once. The built-in derived fields are all of 
type FELjnap.then.inCerpolate.derived-f ield*. This option produces 
the same results one gets with precalculated (eagerly evaluated) component fields. 
Note that for vertex-based queries {i.e., at-cell()), which entail no interpola- 
tion, the distinction between FELjnap.then_interpolate_derived_f ield* 
and FEL.interpolate.thenjnap.derived.f ield* is immaterial. 


19.4 Built-in derived fields 

FEL provides several “prepackaged” derived fields. The differential operator fields are 
specialized derived fields which are discussed separately in Chapter 20. The rest of the 
FEL built-in derived fields are each briefly described in the following list. The sam- 
ple declarations indicate the type of the derived field produced, and also the type(s) of 
the component fields which must be provided as arguments. The generic component 
fields float-field and vector3f. field represent fields already created in an 
application, from which you want to derive some new fields. The newly created de- 
rived fields are called derived-f loat.f ield^nd derived_vector3f -field. 
All derived fields are templatized, but for clarity the typedef’d names are used here. 
Therefore further variations on these built-in derived fields can be generated by using 
the templatized declarations directly. 

• FEL-f loat-f ield.ptr derived-float-field = 

new FEL-magnitude-of.vector3f-f ield (vector3f_f ield) ; 
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derivecLf loat.f ield now points at a newly created scalar field whose val- 
ues are the lengths (Euclidean norms) of the vectors at the corresponding points 
in vector3f.f ield. 

• FEL_float.field.ptr derived_f loat.f ield = 

new FEL.absolute.value.of.float.fieldf float-field) ; 

derivecLf loat.f ield now points at a newly created scalar field whose val- 
ues are the absolute values of the scalar values at the corresponding points in 
f loat .field. 

• FEL_f loat.fi eld.ptr derived-f loat_f ield = 
new FEL_negate_of.float_f ield ( float-field) ; 

derived-float-field now points at a newly created scalar field whose val- 
ues are the additive inverses of the scalar values at the corresponding points in 
float-field. 

• FEL-vector3f.field.ptr derived-vector 3 f_f ield = 
new FEL_negate.of.vector3f.f ield(vector3f.f ield) ; 

derived-vector 3 f .f ield now points at a newly created vector field whose 
values are the additive inverses of the vector values at the corresponding points 
in vector-field- that is, the corresponding vectors in vector.f ield and 
derived_vector_f ield have equal magnitudes but opposite directions. 

• FEL.f loat.f ield-ptr derived-float .field = 

new FEL_suiu-of-f loat.f ield ( f loat.f ieldl , f loat.f ield2 ) ; - 

derived.f loat.f ield now points at a newly created scalar field whose 
values equal the sum f loat.f ieldl and f loat.f ield2 at corre- 
sponding points - that is, derived.f loat.f ield = f loat.fieldl + 
f loat.f ield2. 

• FEL.vector3f.f ield.ptr derived.vector3f .field = 
new FEL_sum_of _vector3 f .field (vector3f.f ieldl , vec- 
tor3f.f ield2 ) ; 

derived.vector3f. field now points at a newly created vector field 
whose values equal the vector sum of vector3f.f ieldl and vec- 
tor] f _f ield2 at corresponding points - that is derived.vector.f ield 
= vector 3 f.f ieldl _+ vector 3 f.field2. 

• FEL.f loat.f ield.ptr derived.f loat.f ield = 

new FEL.dif f erence.of. float-field ( float_f ieldl , 
f loat.f ield2) ; 

derived.f loat.f ield now points at a newly created scalar field whose val- 
ues equal the differences between f loat.f ieldl and f loat.f ield2 at cor- 
responding points - that is, derived.f loat.f ield = f loat.f ieldl - 
f loat.f ield2. 
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Customizing the built-in derived fields 
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field(s), and the TO parameter represents the type that is derived from them. These 
parameters must be replaced with the actual input and output types of the derived field 
you wish to construct - using the appropriate syntax, of course. 

For example, say you want to multiply a vector field by a scalar field, that is, 
you want a derived vector field which scales the vectors of one of its component 
fields by the scalar values of its other component field. Here’s a declaration using 
FEL_product-f ield: 


FEL_vector3f_f ield_ptr scaled_vector_f ield = 
new FEL_product_f ield<FEL_vector3 f , FEL_vector3f , f loat> 
(unscaled_vector_f ield, f loat_field) ; 

The arguments in the <>’s are the template parameters. The first such argument 
is the TO type: this derived field will produce vectors. The second and third template 
parameters are the FROM types - in other words, the types of the component fields from 
which the derived field will produce its values. 

The arguments in the ( ) ’s are fodder for the derived field’s constructor. These 
arguments are pointers to fields of type FROMl and FROM2 — in this case, an 
FEL_vector3 f.f ield.ptr, and an FEL-f loat-f ield_ptr, respectively. These 
fields must both share the same mesh, and must be properly initialized at the time the 
derived field is constructed. All the constructors for the built-in derived fields require 
pointers to fields of the FROM types (and, optionally, a name for the field). 

In this example derived field (scaled-vector-f ield), the vectors of un- 
scaled.vec tor. f ield are scaled by the floating point values of float-field at 
corresponding locations. One could invoke a single (location invariant) scaling factor 
by making float-field an FEL.constant.f ield (see Chapter 17). 


19.5 PLOT3D derived fields 

In addition to the generic built-in derived fields described in the last section, FEL sup- 
plies over fifty specific derived fields defined by PLOT3D. These predefined derived 
fields can be created by calls to specialized convenience functions, or by enum’ed re- 
quests to a PLOT3D “field manager”. For more information on both of these options, 
see Chapter 24. 


19.6 Constructing a custom derived field 

If the built-in derived fields don’t offer the functionality you need, you can define your 
own from scratch. And even if you can achieve the same functionality by other means, 
encapsulating certain data transformations in a derived field can be an effective pro- 
gramming strategy. The most important and involved step in defining a derived field is 
creating the mapping function. 
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19.6,1 Writing a mapping function 

The mapping function defines how derived field data are actually derived from the com- 
ponent field(s). The mapping function must be declared and defined by the user, and 
passed to a derived field constructor as a pointer. Each of the derived field constructors 
expects a pointer to a function with one of three prototypes, depending on whether the 
derived field operates on one, two, or three component fields: 

int map_func (const FEL._solution_globals&, 
const FROM1 * , 
void* , 

TO*) 

int map_func ( const FEL_solution_globals&, 
const FROMl* , const FR0M2*, 
void* , 

TO*) 

int map_func (const FEL_solution_globals&, 

const FROMl*, const FROM2 * , const FROM3 * , 
void* , 

TO*) 

Of course, you can name your own mapping function whatever you’d like; 
“map June" is used here solely as an illustration. 

TO, FROMl , FR0M2 , and FR0M3 are template “type placeholders”, which are 
replaced by actual types in actual mapping function declarations and definitions. The 
only difference among the three prototypes is the number of FROMs, which corresponds 
to the number of component fields in the derived field. In all cases, the TO parameter 
represents the node type of the derived field itself. Thus, the FROMs represent the input 
to the mapping function, the TO represents the output, and the guts of the mapping 
function just specify how the FROMs produce the TO. 

Here is an example of a mapping function, which calculates the magnitude of the 
difference between two input vectors. This can be used to construct a derived scalar 
field, whose isosurfaces, for instance, show the spatial pattern of discrepancy between 
two vector fields. 

int vector_difference_mag (const FEL_solution_globals&, 

const FEL_vector3 f * vO, 
const FEL_vector3f * vl, 
void* , 
float* mag) 

{ 

*mag = FEL_magnitude ( *v0 - *vl) ; 
return 1; 


} 
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This function just subtracts one input vector (vl) from the other (vO) and passes 
the vector result to FELunagnitude, whose floating point return value is stored into 
the location pointed to by mag. The FEL_solution_globals& and void* are 
unnamed here to preclude compiler warnings about unused variables. The first three 
arguments must be declared const to inform the compiler that they are not modi- 
fied by the mapping function. You must be sure to include these const declarations, 
and you must be sure to honor them (i.e., don’t try to modify these arguments in the 
mapping function), or your code will not compile. Also note that the first argument 
is a reference and must be so declared. Omitting either the const qualifier or the 
reference declarator (“&”) will usually cause a compiler error in the derived field con- 
structor taking this mapping function (with a message like “no instance of constructor 
[for whatever derived field you are trying to make] matches the argument list”). 

The transformation of FROMs to TO may involve any sort of mathematical gymnas- 
tics one wishes. Additional data which may (or may not) figure in the transformation 
are available from the FEL-Solution.globals and (if you provide it) the void*. 

Using the FEL-Solution.globals 

The FEL_solut ion.globals contain a collection of values associated with the de- 
rived field. Four of these values - 

float f r ee_s t ream_mach ; 

float alpha; // angle of attack 

float reynolds_number ; 

float time; // usually not used 

— are ultimately taken from the PLOT3D solution file header of the first component 
field (or one of its ancestors, if it is a derived field itself). Since the multiple component 
fields of a derived field must share a common mesh, they will typically also share the 
same header data; but be careful if you are somehow combining disparate data. There 
are a handful of other fields in the FEL.solution.globals structure which are 
used by various built-in mapping functions. 

Here is another sample mapping function, which is set up to receive a veloc- 
ity vector input and uses the the f ree.streanumach and alpha values from the 
FEL-Solution_globals to calculate the perturbation velocity: 

int perturbat ion_velocity ( const FEL__solution__globals& sg, 

const FEL_vector3f * v, 
void*, 

FEL_vector3f * pv) 

{ 

FEL_vector3f vinf; - - 

vinf[0] = sg . f ree_stream_mach * cos ( (M_PI/180 . ) * sg. alpha) 
vinf [1] = sg. f ree_stream_mach * sin ( (M_PI/180 . ) * sg. alpha) 
vinf [2] =0.; //no sideslip angle beta ... 

*pv - *v - vinf; 
return 1? 


} 
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Client data via void* 

The void* can be used to provide arbitrary client data to the mapping function. At 
creation time, you can register a void* with the derived field, and this pointer will 
be passed to the mapping function every time it is invoked. This mechanism provides 
the mapping function with directions to find any additional data it may require - data 
which can be specifically updated for a given query on the derived field. 

An example using client data passed to a mapping function will be given below. 

Handling exceptional cases 

Mapping functions may encounter problematic data. For example, a mapping function 
may receive a “0” input value which figures in the denominator of some expression. 
Singular cases like these can be indicated by the int return value of the mapping 
function. The derived field at_calls forward this return value to the client. As usual, 
the client should check the return value of the at-call, and be prepared to respond 
appropriately. Ignoring the return value is perilous: if the at_call fails, the field 
value “return slot” will not be updated, and therefore will contain stale and possibly 
misleading data. 

19.6.2 Derived field declarations and constructors 

There are two types of derived fields, which allow you to determine the relative or- 
dering of mapping and interpolation (see Section 19.3). If you want mapping to pre- 
cede interpolation, use FELjnap-then_interpolate.derived_f ield*. If you 
want interpolation to precede mapping, use FEL-interpolate-thenjnap-der- 
ived.f ield*. 

The *’s in these class names are wildcards for the number of component fields of 
the derived field. The * should be replaced by “1”, “2”, or “3”, as appropriate. 

The non-built-in derived fields must be declared using template syntax. The tem- 
plate arguments are enclosed in angle brackets (<>s), immediately following the de- 
rived class type specifier. The template arguments specify the type of the derived field 
itself (the output type of the mapping function), and the type(s) of the component 
field(s) (the input type(s) of the mapping function). The derived field type (TO) is 
given first, followed by the component field type(s) (FROM). 

A derived field which implements a “dot product 1 ’ operation takes two vector fields 
(FROM! and FROM2) and uses them to calculate a floating point value (TO). The tem- 
plate part of the declaration would look like: 

FEL_map_then_interpolate_derived_field2<f loat , FEL_vector3f , 

FEL_vector3f > 

or 

FEL_interpolate_then_map_derived_field2<f loat, FEL__vector3f , 

FEL_vec tor 3 f > 
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Following the angle-bracketed template arguments is a parenthesized list of argu- 
ments for the derived field’s constructor. In order, these are a list of pointers to the 
component fields, a pointer to the mapping function, a pointer to client data, and an 
optional name. The component field pointers should be listed in the same order as they 
appear in the mapping function — this ordering establishes the correspondence between 
component field values sent to and received by the mapping function. In both C and 
C++ using to get the address of a function is optional, so the name of the mapping 
function will serve as a pointer; however, by all means use if it makes you feel 
better. 

Assuming that we already have on hand component vector fields in- 
put.vector-f ieldl and input-vector_f ield2, and a mapping function 
make.dot.product that converts two vectors into a float, complete derived field 
declarations follow: 

FEL_f loat_f ield_ptr dot_product_f ield = 

new FELjmap_then_interpolate_derived_f ield2 
<f loat , FEL_vector3f , FEL_vector3f > 

( input_vector_f ieldl , input_vector_f ield2 , 
make_dot .product, NULL) ; 

FEL_f loat_f ield_ptr dot_product_f ield = 

new FEL_interpolate_then_map_derived_f ield2 
<float, FEL_vector3f , FEL_vector3f > 

( input_vector_f ieldl , input_vector_f ield2 , 
make_dot.j?roduct , NULL, "my_dot_product_f ield" ) 

The ' ' NULL ' ' argument fills the slot for the (unused) void* client data. In 
the second example, we’ve given the derived field an (optional) name, which can be 
retrieved via dot.product.f ield->get .name ( ) . It often improves readability 
to break up the template and constructor arguments over several lines, as we’ve done 
here. A healthy compiler won’t mind. 

19.6.3 Derived field checklist 

Declaring and constructing a derived field requires some planning. Here is a little 
checklist: 

• Map then interpolate, or interpolate then map? The relative or- 

dering of these two operations is determined by the type of de- 
rived field you construct; i. e. , there are separate types of derived 
fields (FEL_map-theruinterpolate-derived_f ield* and 

FEL_interpolate-then_map-derived-f ield* for each option. 

See Section 19.3 above. 

• How many component fields are there? Component fields are the “input” to the 
derived field, from which its values are derived. FEL allows one, two, or three 
component fields. This number immediately follows “field” in the complete 
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derived field class name. It is unusual to need more than three component fields. 
If you have more than this number of component fields, you may be able to 
“factor” the derivation into parts, each involving three or fewer components, and 
construct the overall derivation incrementally, 

• What are the component field types? These types (typically float, 
FEL_.vector3f, or FEL_plot3cLq) are needed for the template arguments 
in the derived class declaration. 

• What are the component field instances? The component fields (of the types just 
mentioned in the item above) on which the derived field is built need to exist 
when the derived field is constructed. FEL pointers to these fields are passed as 
arguments to the derived field constructor. 

• What is the mapping function? The mapping function converts component field 
types to the type of the derived field itself. A pointer to this function must be 
provided to the derived field constructor. (Or at least one must provide a func- 
tion pointer of the appropriate type, which should actually point to a function of 
the approriate type when the derived field is queried. By redirecting the func- 
tion pointer registered with the derived field, one could dynamically change the 
derivation, based on some conditions, just prior to querying the field...) 

• Are there any client data which need to be passed to the mapping function? 
A pointer can be registered with the derived field at construction time, which 
provides arbitrary data to the mapping function. The data can be changed at any 
time, for instance, just prior to querying the field. 

19.6.4 A more or less complete derived field example 

Here is an example which constructs a derived field representing the signed distance 
from a plane. Such a field could be used in conjunction with an iosurface routine, for 
instance, to implement a cutting plane routine. 

typedef struct 

{ 

FEL_vector3f n; // normal 

float p; // point 

} hessian; // Hessian normal form of a plane 

// the mapping function 

int signed_distance (const FEL_solution_globals5c, 

const FEL_vector3f& test_point, 
void* plane, // client data 
float* dist) 

{ 

hessian* hnf = (hessian* ) plane; 

// calculate signed distance 
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*dist = FEL_dot (hnf->n, test_point) + hnf->p; 
return 1; 


int cutting_plane (FEL_plot3d_f ield__ptr plot3d_f ield, 

hessian* plane, 

FEL_f loat_f ield_enum scalar_f ield) 

{ 


FEL_mesh_ptr mesh = plot3d_field->get_mesh() ; 

// "convert" mesh to position vector field 

FEL_vector3f_f ield_ptr pvec = 

new FEL_mesh_as_vector3f_f ield (mesh) ; 

FEL_f loat__f ield_ptr distanced ield = 

new FEL_map_then_interpolate_derived_f ieldl 
< float , FEL_vector3f > 

(pvec, signed_di stance, plane); 

FEL_f loat_f ield_ptr color_by_f ield = 
plot3d_field->make_float_f ield (scalar_f ield) ; 

// colormapping via "color_by_f ield" 
isosurf (distance_f ield, color_by_f ield, . . .) ; 

} 

In this example, the typedef d struct hessian contains a point and normal defin- 
ing a plane. The mapping function signed-distance { ) uses the point and normal 
to calculate the signed distance from the plane (negative on one side, positive on the 
other) to an incoming test .point. The mapping function is used by the derived 
field distance-field, along with the plane equation, which is provided as client 
data (actual argument plane) to the derived field constructor. The derived field is built 
on top of a single component field — pvec — which is a field of the position vectors of 
the underlying mesh, courtesy of the adaptor class FEL_mesh.as.vector3 f .field 
(see Chapter 17). When distance-field is queried for a value at a given location, 
it returns a signed scalar indicating the distance of the query location from the plane 
described in plane. Thus an at.cell ( ) call on distance-field will return an 
array of values indicating the distance of each vertex of a cell from a plane. Interpo- 
lating the 0-valued surface on this field (via marching cubes, or the like) will yield the 
cutting plane specified by plane. 
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Figure 19.1; A visualization illustrating several uses of derived fields. The white line 
originating from the trailing edge of the shuttle wing is a streamline traced in the ve- 
Tocit\ derived field. (Mng the streamline as a local frame of reference, cutting planes 
were a nstructed as in the example in this chapter, and then each plane was trimmed 
to a d:>k. Textured on each disk is the pressure derived field. Visualization by Chris 
Henze 



Chapter 20 

Differential Operator Fields 


20.1 Gradient, divergence, and curl 

The FEL differential operator fields are specialized derived fields which apply the dif- 
ferential operators grad, div, and curl, to a preexisting base field, in three-dimensional 
Euclidean space. The three operators are typically written in terms of a single operator, 
the nabla operator (V). 


20.1.1 Grad 


The gradient operator acts on a scalar field and produces a vector field indicating the 
local direction and magnitude in which the scalar field is changing most rapidly. In 
Cartesian coordinates, where / is a scalar field, f = f[x, y, z ): 


ffra d/ = V/=gi + gj + |{k 


The gradient vectors are always perpendicular to the level surfaces (isosurfaces) of /. 
The spatial derivative of / in any direction is just the projection of the gradient vector 
on that direction. 


204.2 Div 


The divergence operator acts on a vector field and produces a scalar field indicating the 
local net outward flow per unit of time and volume. In Cartesian coordinates, where F 
is a vector field, F = F(x, y, z) r and F x is the x-component of F (likewise for y and 

zy 


„ ^ dF x dF % 
div F = V- F = ir ^ + - 7r 1 
dx dy 


dz 


Locations in a vector field F where V • F > 0 are called sources (more outflow than 
inflow), locations where V • F < 0 are called sinks (more inflow than outflow), and 
locations where V * F = 0 are called source-free. If the entire field is source-free, F is 
called solenoidal. 
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jFEl_aradlenLrwld1<T0,FR0M> | 
j |FEL_differentiaropo<'atorJield1<TO.FROM>i ^— |FEl_dfvergence_field1<T0,FR0M>] 


|FEL_diffefentlai_operatorJlekl<TO,FROM>| ( 


lFEL_curi_fteld1 <TO,FROM> 1 
> jFEL _flradient_fMd2<TO t FROM>l 
V iFEl_differential_operator_field2<TO,FRQM> ^HFEL_divergence_fletd2<TO,FROM>l 


lFEL_curi_ffeid2<T0,FR0M> j 


Figure 20. 1 : The FEL differential operator field class hierarchy. 


20.1.3 Curl 


The curl operator acts on a vector field and produces another vector field indicating the 
local direction and magnitude of the rotation of the original vector field. In Cartesian 
coordinates, where F is a vector field, F = F(x, y, z ), and F x is the x-component of 
F (likewise for y and z ): 


curl F = V x F = 


(0J\ 

V dv 




k 


At a given location in a swirling vector field, the curl is proportional to the local angular 
velocity. The curl is nonzero in “straight” vector fields with shear. Vector fields with 
V x F = 0 everywhere are called irrotationaL 


20.2 First-order and second-order accuracy 

As you can see from the descriptions in the previous section, all FEL differential op- 
erator fields require spatial derivatives of their scalar or vector base field values. In 
FEL, these spatial derivatives are estimated by either a first-order accurate scheme, or a 
second-order accurate scheme, depending on the type of differential operator field you 
construct. Figure 20. 1 shows that the FEL differential operator field class hierarchy has 
two parallel lineages: the “1” or “2” in the class names determines whether a first-order 
or a second-order accurate scheme is used in estimating the derivatives used by a given 
differential operator..;: - — ,, : . — .•*, ' 1 - 

In the first-order differential operator fields (descendants of FEL_differ- 
ent ial.operator.f ieldl<TO, FROM>), the same interpolating polynomial 
which is used to interpolate values on the base field is analytically differentiated to pro- 
duce an expression yielding interpolated derivatives. If the differential operator field in- 
terpolation mode (see Chapter 9) is FEL-ISOPARAMETRICJCNTERPOLATION, the 
computational space shape functions on a given cell are differentiated with respect to 
the computational coordinates, the resulting partial derivatives are evaluated at the ap- 
propriate locations, and then transformed into physical space by way of the metrics. If 
the interpolation mode is FEL.PHYSICAL_SPACE.INTERPOLATION, the physical 
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space interpolating polynomial on evatoted directly 

ica, space * Ve resulting physical space partial d '"” 1 ' veS are “ 

Mds (descendants of 
In the second-order dtfTerennd pe FR01)>)i partial den»attves 
FEL.dif Eerer.cia 1 -Operator differencing in computational space, 

are estimated a, each vertex ^ b, dt. medics. Then die 

and these derivatives are "“**“1““^” e ,ch cell vertex and, if necessary, 
differential operator quantities ar g 

interpolated at interior lo " ati ^, operator fields will produce more accurate 

In general, the second-order different P secon d-order fields require more 

and reliable results than the first-order fie . * field values from a large r stencil, 

work, since the central difference supported on unstruc- 

At present, second-order differentia P ^ ^ differcnce scheme on 

20.3 Creating differential operator fields 

Differential operator fields node^pf^t^ti^t ^veedte^d 

as base fields. As long as the has. t fied to . created in a „y number 

"tf SS.'SK? * »'« “• a derived “• * ^ 

diff S« few L pfided to the di ^“ al d X r e« !a“ Wd°cons« 
pointer. In fact, the only other ar 6“ me string . You cho ose between first- and 

is an optional name, in the orm o a iadn a differential operator field with a 
second-order numerical schemes y haracter in the class name. 

(first-order) or “2” (second-order) as the ^ (parameterize d) by two types: the 

,ypetefrr^S-. T of.^ 

Ma ulf ,he pro “ ,pi “' d,ff '""“ a ' opera ' 

tor field declaration looks like this. 

FEL_pointer< FEL typed field<TO>^> gr.^fxeXd - 

nd FEL^ointer< n rEL^typed_field<FROM> >1 , 
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sr a “ iKeif - remmed * “»• “ i— »y . 

Thus, 

FEL_float_field_ptr temperature_f ield; 


FEL_vector3f_f ield_ptr grad_of_temperature field - 
( temperature_f ield) ; 


t ureTf i e 1 c^Usi n^a^tynedef th *“* retUmS ** gradient of tempera- 

• in 8 a typedef, the same declaration could be accomplished as; 

Ptr tempera ture_f ield; 




r “ le ^ c *— ptr grad_of_temperature field - 
FEL - gradient _of_float_f ield (temperature^.f ield) , 

Here are a few more annotated examples. 

FEL_ flo at_ f ieid_ptr temperature_field; 
EL - vector3f - f ield velocity_f ield; 


/ / isc order 


FEL_float_field_j?tr div_6f__velocity field = 

new FEL_divergence_ 0 f_vector3f_f ieldltvelocity_f ield) ; 

/ / 2nd order 

FEL_float_field_ptr div_of_velocity field = 
new FEL_diver g ence_of_vector3f_f ieTd 2 tve l0 ci t y_ f iel dj ; 

pp T 2 JJ d °f der ' double precision output 
FEL_double_field_ptr div_of_velocity_f ield = 
new FEL_divergence of vector3f kl 

{velocity_f ield) 7 ector3f - f ield2<double, FEL_vector3f > 

II curl of velocity = vorticlty ~ : 

FEL_vector3f_f ield_ptr vorticity.f ield - 
new FEL_curl_f ield2<FEL_vector3f , FEL_vector 3 f> (velocity.f ie 

II typedef 'd version 

FEL_vector3f_fieId_ptr vorticity_f ield = 

new FEL_curl_of_vector3f_field2(velocity_field) ; 
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20.4 “Chaining” differential operator fields 


As mentioned previously, the base field of a differential operator field can itself be a 
differential operator field. This allows operator “chaining” to create new operators. For 
example, the Laplacian operator (V 2 , sometimes written A), in Cartesian coordinates, 
operating on a scalar field / = /(x, y, z), is defined to be: 


A / = div grad f = V 


~ dx 2 + dy 2 + dz 2 


In light of this definition, one can declare a scalar Laplacian differential operator field 
as follows: 


FEL_float_.fi eld _ptr pressure; // a derived field: de- 
fined elsewhere . . . 


FEL_vector3 f_f ield_ptr grad_pressure = 
new FEL_gradient_of_vector3f_f ield2 (pressure) ; 

FEL_f loat_f ield_ptr laplacian_pressure = 
new FEL_divergence_of_vector3 f_f ield2 (grad_pressure) ; 

Now laplacian.pressure is a scalar field which returns the second spatial 
derivative, or “curvature”, of the pressure field. Note that pressure is itself a derived 
field, whose values are only calculated on demand. 

Chaining the differential operator fields is a powerful technique, and FEL supports 
chains of arbitrary length; but one must be aware of the numerical limitations of grid- 
ded data. Numerical differentiation by way of finite difference schemes magnifies noise 
in the data, and is subject to truncation and roundoff error. These factors can quickly 
overwhelm any meaningful signature in a dataset, due to the repeated numerical differ- 
entiation entailed by a chain of differential operator fields. 

The numerical problems associated with repeated differentiation are particularly 
severe on unstructured grids. Unstructured grids support only first-order differential 
operators, and these produce constant values across a given tetrahedron. Subsequent 
differentiation at the nodes of the constant-valued tetrahedra invariably yields zero 
derivatives. On structured grids, chained differential operator fields should be second- 
order, to avoid the derivatives from “bottoming out” prematurely. 

In addition to these numerical issues, there are performance considerations asso- 
ciated with chained differential operator fields. In the second-order central difference 
scheme, evaluating the derivatives at a given vertex requires base field values from a 
neighborhood of adjacent vertices. If the base field itself requires central difference 
values, it will have to fetch values from yet a wider neighborhood, encompassing the 
first. This “expanding neighborhood” (or “stencil”) around each vertex results in signif- 
icant overhead for a field query. In addition, adjacent vertices have largely overlapping 
“expanded neighborhoods”, so adjacent field queries result in a lot of duplicate work. 
For these reasons, if performance is an issue, it may make sense to eagerly evaluate 
(see Chapter 17) or at least cache (see Chapter 17) the results of a chained differential 
operator field, particularly if it involves highly derived fields. 
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Here is an example combining differential operator fields and derived fields. 
The Laplacian operator, in Cartesian coordinates, operating on a vector field F = 
F(x,y, z ), is defined to be: 

d 2 F d 2 F 

AF = grad div F - curl curl F = V(V-F)-Vx(VxF)= + 


We can define a vector Laplacian operator field as follows: 

FEL_vector3f_f ield_ptr v_field; // some vector field, de- 
fined elsewhere 

FEL_f loat_f ield_ptr div_field = 
new FEL_divergence_of_vector3f_f ield2 (v_f ield) ; 

FEL_vector3f_f ield_ptr grad_div_f ield = 
new FEL_gradient_of_float_f ield2 (div_field) ; 

FEL_vector3f_f ield _ptr curl_field = 
new FEL_curl_of_vector3f_f ield2 (v_f ield) ; 

FEL_vector3f_f ield_ptr curl_curl_f ield = 
new FEL_curl_of_vector3f_f ield2 (curl_f ield) ; 

FEL_vector3f_f ield_ptr vector_laplacian_f ield = 
n ew FEL_di f f e r enc e_o f _vec t or 3 f _f i e Id 
(grad_div_f ield, cur l_curl_f ield) ; 

// if desired, convert into core field 

vector_laplacian_f ield = 
vector_laplacian_f ield->get_eager_f ield( ) ; 



Chapter 21 

Instantiating Fields 


The most prominent use of templates in FEL is for the node type of fields. The tern- 
plating makes it possible for the user to construct fields with a new node type with a 
minimal amount of code. FEL requires a few basic operations be supported for the 
node type, so that the library can interpolate if necessary when queried about field val- 
ues. The number of operators required is intentionally kept small to make it easier to 
introduce custom types. 

To demonstrate the minimal requirements of the node type, we present a few ex- 
ample types below. Keep in mind that most built-in numerical types support all the 
required operations, and then some. 


21.1 Basic type requirements 

The first example type is the “f oo” type. A f oo has basically the behavior of a spartan 
scalar. 

class foo { 
float value; 
public; 

food { } 

foo(float v) : value (v) { } 

friend foo operator* (double d, const foot f) { 
return foo( (float) (d * f.value)); 

} 

friend foo operator-*- (const foo& lhs, const foo& rhs) { 
return foo ( lhs .value + rhs. value) ; 

} 

// ostream operator not necessary, but handy 
friend ostream& opera tor« (os t reams strm, const foo& f) { 
return strm << f.value; 

} 

} ; 
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A “f oo” must have a default constructor (so that one can allocate an anay of them), 
and we provide a constructor with a float argument so the example is not too trivial. 

The * and + operators support multiplying a f oo by a scalar and adding two f oo 
objects together. In a program, the instantiation of a f oo field would look like: 

main() { 

/ / make a dummy mesh 

FEL_mesh_ptr mesh = new FEL_regular_mesh ( 5 , 7, 11); 

// convenience typedefs for foo fields 

typedef FEL_f ield<f oo> FEL_f oo_f ield; 

typedef FEL _jpointer<FEL_foo_f ield> FEL_f oo_f ield_ptr ; 

typedef FEL_core_f ield<f oo> FEL_core_foo_f ield; 

// make a foo field and iterate over it 
foo* foo_data = new foo [mesh ->card (0) ] ; 

// should fill in the foo data buffer ... 

FEL_f oo_f ield_ptr foo_field; 

foo_field = new FEL_core_f oo_f ield (mesh, foo_data) ; 
FEL_vertex_cell_iter iter; 
int res; 

for (foo__f ield->begin (&iter) ; ! iter . done () ; ++iter) { 
foo f; 

res * foo_f ield->at_vertex_cell ( *iter , &f ) ; 
assert (res == FEL_0K) ; 

cout << "field value at " << *iter << " is " << f << endl; 

} 

The typedef statements are not essential, but they come in handy further down the 
line, since the template syntax can get a bit tedious. 

An example closer to what one might do in practice involves the construction of a 
field where each node has a vector of 10 doubles: 

typedef FEL_vector<10 , double> FEL_vectorlOd; 
typedef FEL_f ield<FEL__vectorlOd> FEL_vectorlOd__f ield; 
typedef FEL_pointer<FEL_vectorlOd_f ield> FEL_vectorlOd_f ield_ptr 
typedef FEL_core_f ield<FEL_vectorlOd> FEL_core_vectorlOd_f ield; 
FEL_vectorlOd* vectorlOd_data = 
new FEL_vectorl0d[mesh->card(0) ] ; 

FEL_vectorlOd_f ield_ptr vectorlOd_f ield = 

new FEL_core_yectorlOd_ field (mesh, vectorlOd_data) ; 
for (vectorlOd_f ield->begin (&iter) ; I iter . done ( ) ; ++iter) { 
FEL_vectorlOd v; 

res = vectorlOd„f ield->at_vertex_cell ( *iter # &v) ; 
assert (res == FEL_0K) ; 

cout << "field value at * « *iter « " is " « v « endl; 

} 
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} 

The vector template is defined by FEL (see Chapter 3); there is no need to provide the 
required math operators since they are already in the library. 


21.2 Differential operator field requirements 

If differential operator fields will be constructed with the new node type, then a few 
more operators must be defined: 

class bar { 
float value; 
public : 
bar { ) { } 

bar (float v) : value (v) { } 

friend bar operator* (double d, const bar& b) { 
return bar ((float) (d * b. value) ); 

} 

friend bar operator* (const bar& lhs, const bar& rhs) { 
return bar ( lhs .value + rhs. value); 

} 

friend bar operator- (const bar& lhs, const bar& rhs) { 
return bar (lhs .value - rhs. value); 

} 

// the following should not be necessary, 

// but are needed for some SGI compilers, which can 
// get confused at instantiation time 
friend bar operator* (const bar& b, double d) { 
return bar(b. value * d) ? 

} 

friend bar operator- (const bar& b) { 
return bar ( -b . value) ; 

} 

friend bool operator== (const bar& lhs, const bar& rhs) { 
return lhs. value == rhs. value; 

} 

} ; 

Technically speaking, the only extra operator that should be required is bar 
operator- (const bar&, const bar&) , Unfortunately, as noted in the ex- 
ample, a few more operator definitions may be necessary if the compiler and linker 
that one is working with get confused. Fortunately, the error messages indicating the 
need for a particular operator are typically not too difficult to decipher* and the opera- 
tors are fairly straight-forward to define. 
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Currently FEL contains typedef statements for fields with the following node 
types: 

• float 

• double 

• FEL_vector2f 

• FEL_vector3f 

• FEL_vector3d 

• FEL-matrix33f 

• FEL.plot3d-density-momentum 

• FEL_plot3d_q 

The last two types are aggregates defined for the solution vector data standard to 
PLOT3D [WBPE92]. The “q” type contains the entire PLOT3D solution, the “den- 
sity-momentum” field contains the subset of the variables necessary for velocity- 
related derived fields. The math operators for both types simply define the same oper- 
ators in turn for each component in the object. So, for instance, operator + with two 
q objects returns a new object where the density, momentum, and energy components 
are each the sum of the corresponding components in the two arguments. 


Chapter 22 

File I/O 


To construct instances of most mesh and field classes in FEL, one provides, as argu- 
ments, parameters and pointers to buffers containing data. Typical parameters include 
specification of structured mesh dimensions, and data buffers typically include vertex 
coordinates or solution data. For most scientists, mesh data and solution data are stored 
in files of various formats. FEL provides file reader functions that extract the data from 
files and, where appropriate, that load the data into main memory with a layout ap- 
propriate for the class to be constructed. Using the reader functions, one can simply 
specify the file name and optional flags describing the file format, and the appropriate 
objects will be constructed. 

FEL has families of file reader functions for several major file formats. The most 
extensive support in FEL is for the PLOT3D [WBPE92] format. A second family of 
readers accepts paged PLOT3D files (see Chapter 23). Paged PLOT3D files contain 
the data of standard PLOT3D files, but the data are reorganized into page-sized chunks 
that can be read in on demand. FEL contains less extensive support for two more file 
formats: FITS [FIT] and Vis5D [Vis]. Finally, the library provides a generic set of 
reading functions that are independent of a particular file format. The idea is that an 
application written in terms of the generic routines can work with a variety of file for- 
mats. The generic routines currently support work with PLOT3D and paged PLOT3D 
files, with limited support for the Vis5D format. 

File readers are built using standard FEL operations. Thus, it is possible to write a 
reader for a new file type without having to modify FEL. 


22.1 PLOT3D file reader functions 

The visualization application PLOT3D [WBPE92] defines a family of file formats for 
storing meshes and fields. The FEL file reader functions handle files designated as * 3D 
(/WHOLE)” by PLOT3D. The PLOT3D “ID”, “2D”, and “3D (/PLANES)” formats 
are not supported. This section contains descriptions for the three PLOT3D-related file 
reader families. Most of this section describes the family of generic file readers, with 
notes describing how the other two families differ. A later subsection shows how the 
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functions correspond among the three families. 

Most applications should use the generic functions so that they work with both 
types of files. While the generic functions have some amount of overhead over the 
specific functions, the overhead is quite minimal and is outweighed by allowing your 
program to work with multiple file formats. Future versions of FEL may enhance the 
generic file reader functions which would let your application read those files with little 
or no effort. 

22.1.1 The PLOT3D flags 

The FEL PLOT3D flags allow the user to specify the type of a file. They spec- 
ify the type of the file (PLOT3D or paged file) as well as the particular variation 
of the PLOT3D file format. The generic file readers use the file type flags (ei- 
ther FEL-PLOT3D_3D_WHOLE or FEL-PAGED.PL0T3 D_3 D.WHOLE) to determine 
the file-type-specific routine that should be called to do the actual work. The variation 
flags tell the standard PLOT3D routines which PLOT3D variation the file contains. 
The paged PLOT3D routines ignore the variation flags because the paged files are self- 
identifying. If you call a standard-PLOT3D- or paged-PLOT3D-specific reader routine 
directly, the flags argument must correspond to the file type that the routine handles. 

Several flags can be expressed simultaneously using the C 44 1” bitwise-or op- 
erator. For example, “FEL-PLOT3 D_3 D WHOLE | FEL_PLOT3D_MULTIZONE j 
FEL.PLOT3 D_I BLANK” specifies a PLOT3D multi-zone file which includes IBLANK 
information. 

The PLOT3D readers recognize the following flags: 

FEL-PLOT3D.3 D.WHOLE for all PLOT3D (non-paged) files 
FEL-PAGED-PL0T3 D_3 D.WHOLE for PLOT3D paged files 
FEL-PL0T3D.IBLANK file contains IBLANK information 
FEL.PL0T3 D-MULTI ZONE file has multiple zones .... ... . 

FEL-PLOT 3D -UNSTRUCTURED file contains an unstructured (tetrahedral) mesh 
FEL_PLOT3D_FORTRAN -UNFORMATTED binary with FORTRAN control words 
FEL-PLOT3D-FORTRAN -FORMATTED PLOT3D ASCII format 
FEL-PL0T3 D-FUNCTION file is a function file 
FEL-PL0T3DXITTLE.ENDIAN bytes are in little endian order 

22.1.2 Automatic mesh type deduction 

FEL provides an automatic format deducing function, FEL.deduce jnesh.type { ) , 
which makes it easier for the user to read meshes without having to remember the 
particular variation of the file format at hand. The deducer will determine the type of 
the file (standard PLOT3D or paged), and the file’s PLOT3D variation. The deducer 
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allows one to pretend that all PLOT3D files are self-describing. The deducer takes a 
character string file name and returns an unsigned integer representing the flags, or 0 
if it could not successfully deduce the mesh type. The deducer essentially works by 
hypothesizing that the file is a particular format and examining the file to see whether 
the file is consistent with that format. With paged PLOT3D files, the deducer just 
examines the first few bytes of the file to see if the file has a valid paged PLOT3D file 
header. 

Deducing a standard PLOT3D file takes more work. The PLOT3D deducer must 
hypothesize a PLOT3D file variation, interpret the first words in the file as if in the 
hypothetical variation, and then compute the size of the file in bytes implied by the 
header words. If the number of bytes implied matches the actual number of bytes in the 
file, then the deducer returns the specific variation. Otherwise, the deducer continues 
to the next variation hypothesis. The technique used by the deducer is not fool-proof. 
For example, if for some reason the file has extra bytes appended to it, then the implied 
count will not match the actual, and the deducer will fail. 


22.1.3 Reading mesh files 

The function FEL_read_mesh ( ) is the function for reading meshes. The function 
takes one required argument, the character string mesh file name, and returns an 
FEL_mesh subclass. If the reader encounters an error, then it returns NULL. Imme- 
diately following the file name argument is an optional argument with the format flags 
for the given file. By default the reader calls the deducer function to compute the flags. 
The third argument allows reading a specific zone from a multi-zone data set. The 
default value of -1 specifies that every zone should be read in; other values specify a 
particular zone to be read. The function’s declaration is: 

FEL_mesh_ptr FEL_read_mesh (char* , unsigned = 0, 

int = -1) ; 

22.1.4 Getting information about structured mesh files 

Three functions read part of a structured mesh file and return information about it. 
These functions abort if they are called on an unstructured mesh file. All three functions 
return a status code of FEL_OK on success and FEL-FAILED on failure. The first 
argument for the functions specifies the name of the file. The second argument specifies 
the type of the file and must be non-zero (use the deducer if you don’t know the type). 

The same information can be retrieved from a mesh once you have read it, using 
get _n-Zones and get-s true tured_dimens ions. If you do not need the infor- 
mation before reading the mesh, it is more efficient to read the entire mesh and then 
retrieve the information from it. 

The functions are: 

int FEL_read_mesh_dimensions (char* , unsigned, int*, 
FEL_vector3i**) ; 

int FEL_read_mesh_n_zones (char* , unsigned, int*); 
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int FEL_read_mesh_zone_dimensions (char* , unsigned, int, 
FEL_vector3i*) ; 

The first function, FEL_read_mesh-dimensions, returns the number of zones 
and the dimensions of each zone. The dimensions are returned in a dynamically- 
allocated array. This array must be deallocated using delete [ ] after you have fin- 
ished using it. The array is returned by setting the FEL_vector3i* that is pointed to 
by the third argument. 

The second function, FEL_read_mesh_n-Zones, returns the number of zones 
by setting the integer pointed to by the third argument. The last function, 
FEL_read_mesh-zone-dimensions, returns the dimensions of the zone specified 
by the third argument by modifying the FEL_vector3i pointed to by the last argu- 
ment. These last two functions call FEL-read_mesh-dimensions and return part 
of what it returns, so it is more efficient to call that function if you need information 
about all the zones, 

22.1.5 Reading solution files 

FEL provides a set of solution reader functions, all with the prefix FEL-read-, 
e.g., FEL_read_density. Each solution reader function take, as a first argu- 
ment, the mesh upon which the solution is based and a character string file name 
as the second argument. The third required argument is the set of flags spec- 
ifying the specific file format. The flags are the same as those used for the 
mesh reading functions, though some flags (specifically FEL-PLOT3D.I BLANK and 
FEL.PLOT3 D .UNSTRUCTURED) are irrelevant to solution files and are ignored. If 
the flags include F EL J? AGED-PLOT 3 D_3D_WHOLE, i.e., if one will be reading from a 
paged file, then the remaining flags are not necessary, as paged files are self describing. 
The data in the paged file overrides options spelled out by the flags. Immediately fol- 
lowing the three required arguments to the solution reader routines there is an optional 
argument specifying a specific zone to read. By default the data for all the zones of a 
multi-zone file are read. 

Unlike the mesh case, the library does not provide a format deducing function for 
solution files. Typically, one can use the flags returned by FEL-deduce_mesh_type 
as an argument to both the mesh and solution reading functions, since both are usually 
the same PLOT3D variation, e.g., “FORTRAN unformatted”. Thus the lack of a so- 
lution format deducer is typically not an inconvenience. In the cases where the mesh 
and solution are not in the same format, one can still use the mesh deducer result for 
reading the mesh, but the flags for the solution reader must be determined by the user 
and provided manually. 

The typical usage of the deducer and reader functions is: 

unsigned flags; 

FEL_mesh_ptr mesh; 

FEL_f loat_f ield_ptr density_f ield; 

flags = FEL_deduce_mesh_type (mesh_f ile_name) ; 
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22.3 The FITS file reader 

FEL provides a simple reader for FITS files. The reader handles data for 2- or 3- 
dimensional regular meshes, with a single float value at each node. The reader takes as 
its first argument the character string specifying the file name. The next two arguments 
each are a pointer to an integer. FITS file headers consist of a sequence of tags, each 
with a corresponding value. Some users have extended their use of the FITS format 
by defining new, non-standard tags. In particular, the tags “XPOS” or “YPOS”, each 
accompanied by an integer, are used by McDonnell-Douglas to signify the lower comer 
of a subimage (where images are represented by a 2-D mesh). The FEL FITS file reader 
function recognizes the “XPOS” and “YPOS” tags, returning the corresponding values 
in the integers pointed to by the second and third arguments of the reader call. By 
default “XPOS” and “YPOS” are set to 0. 

In general, the FITS reader function does not write any output to the terminal. 
In some cases, it may be handy to request more output, for example, when trying to 
determine why the reader cannot successfully read a particular file. The FITS reader 
routine checks a public global variable, FEL_f it s_reader_verbose, which can be 
set to true by the user. When this flag is true, the reader provides more information 
about what it is reading, including displaying the tags that it encounters in the file. 

22.4 The Vis5D file reader 

FEL provides a simple reader for files in the VIS5D format. There are two routines 
available to the user: FEL_vis5d_get_scalar and FEL.visScLget.vector. 
The scalar routine returns a float field, the vector routine returns a FEL_vector3f 
field. Both routines take a mesh as a first argument and a character string file name as 
a second argument. The scalar field reader takes a third argument: a character string 
providing the name of the variable to read. Finally, both routines take an integer final 
argument specifying the time step to retrieve when working with time series data. 


Chapter 23 


Paged Meshes and Fields 


23.1 Introduction 

Paged meshes and fields are similar to standard FEL meshes and core fields, but they 
use much less memory. Instead of keeping all of the data in memory, they bring in the 
portions that are used. Data are placed into a pool of memory ; the size of the pool can 
be specified by the user. When the pool is full, some data that have not been recently 
used are replaced with the new data. 

A paged mesh or field can be much faster than the corresponding in-memory ver- 
sion if your system does not have enough main memory to hold all of the data. While 
you can rely on the operating system’s virtual memory to bring in the data when you 
are using in-core meshes and fields, the paged versions use memory more efficiently 
and may be able to keep what is currently needed in memory while the operating sys- 
tem cannot. A program will run much faster if its data can be kept in memory. Even if 
you have enough memory, a paged mesh or field can be faster if you only need to use a 
fraction of the data since only that fraction is read from disk. However, a paged mesh 
or field can be slower if you have sufficient memory and repeatedly use all or nearly 
all of the data. This happens because the data are read in smaller amounts compared to 
when the entire file is read at once, and because, for the same amount of data, reading 
data in smaller amounts takes longer than reading it all at once. Also, accessing the 
data is a bit slower with the paged routines because they do some additional checks 
and calculations. More information about the advantages and disadvantages of paged 
files and out-of-core visualization in general can be found in [CE97]. 

From a coding standpoint, the use of a paged mesh or field is very similar to using 
the corresponding in-core objects. The same operations are supported, except that the 
names of some setup calls are different. 

As of this writing, only PLOT3D structured grid and solution files can be paged. 
Support for function files and unstructured grids may be available in a future release. 
Also, the current paging code is not thread-safe, so paged files should not be used in 
multithreaded programs. This should be corrected in the next release. 
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23.2 How paged files work 


The technique that paged files use to bring data from disk as it is used is similar to 
how virtual memory is implemented in computer systems, but using software instead 
of hardware. Like virtual memory systems, the paging code reads fixed-sized blocks 
of data, called pages of data. 

One reason that paged meshes and fields use memory more efficiently than the 
operating system’s virtual memory system is that they can select a different page size. 
The paging code uses pages of 2 kilobytes, while current workstations use page sizes 
of 4, 8, or 16 kilobytes. Smaller pages use memory more efficiently than larger pages 
because there are smaller amounts of unreferenced data surrounding the referenced 
data that must still be read into memory. 

The second reason for paged meshes’ and fields’ more efficient memory usage is 
that they reorganize the data in the file, which places different data on each page. For 
structured grids, an 8x8x8 cube of data from a single zone is placed in a page. The 
original PLOT3D files would place a plane of data on a page. Again, "the reason that 
this file organization uses memory more efficiently is that placing cubes of data on each 
page has smaller amounts of unreferenced data surrounding the referenced data. For 
most 3D traversals of the data, placing 8x8x8 cubes around the referenced data covers 
a smaller portion of the overall data compared to placing (for example) 32x64 planes of 
data around the referenced data. In practice, this reorganization cuts the memory usage 
in half compared to PLOT3D when computing streamlines. If memory is tight, using 
half the memory can translate into larger factors of speed improvement. One drawback 
from reorganizing the data is that the files must be converted to a new format. The next 
section describes how to convert files. 

When you start to use a file, it is opened, and data structures are set up that have an 
entry for each page in the file; each entry indicates that the page has not been loaded 
into memory. When you access some data (e.g., by using an at.cell call), the pages 
where the data are found are computed and checked to see if they are present. If they 
are not present (which would be true if the file was just opened), the pages are read. 
Then, the values are retrieved from the pages and returned. 

Pages read from disk are placed in memory allocated from a pool of pages. There 
is a single memory pool, which means that it is shared among all of the open paged 
files. If all pages are in use when a new page is needed, a page that has not been 
recently used is reused. The size of the pool is user-configurable (see Section 23.5 
below). Preliminary estimates indicate that you can visualize a dataset interactively 
when the pool will hold 20% of the dataset’s grid files and 5% of the solution files. 
These numbers depend on how the data are used and thus will vary considerably. 

There is some memory overhead associated with paged files. When you start using 
a file, the paging code creates data structures that use memory equal to about 0.3% 
of the file size. This can be significant: a 10 gigabyte dataset would use about 30 
megabytes of memory. This memory is allocated separately from the memory pool 
used to hold file data. 
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23.3 Converting PLOT3D files to paged files 

As mentioned above, you need to convert your PLOT3D files to the paged file format. 
The p3dtopaged program will convert one file to the new format. This program can 
be found in the FEL source directory. 

For most files, the program can deduce the type of input file. In that case, the usage 
is: 


p3dtopaged input-file output-file 

The program cannot figure out some file types, especially FORTRAN unformatted 
files with more than 5 solution variables. With these files the type must be specified on 
the command line using flags. The flags are: 

-g indicates that the input is a grid file with no IBLANKS 

-i indicates that the input is a grid file with IBLANKS 

-s indicates that the input is a solution file 

-1 (a one) indicates that the input file has one zone 

-m indicates that the input file has multiple zones 

-b indicates that the input is a binary file 

-f indicates that the input is a FORTRAN unformatted file 

If the program cannot deduce the type of the file, you will have to specify three of 
the above options: one from g, i,or s; 1 orm; and b or f . For example, p3dtopaged 
-imb would specify that the input file is a multi-zone, binary, PLOT3D grid file that 
includes IBLANKS. 

Paged files can be converted to the PLOT3D format with the pagedtop3d pro- 
gram. It only outputs binary files; FORTRAN unformatted is not supported. If your 
original input file was a FORTRAN unformatted solution file with seven parameters 
per node, when you convert the paged file to a PLOT3D file only the first five param- 
eters will be in the file. This happens because paged files contain only the first five 
parameters. If the original PLOT3D file was a binary file, converting the file to a paged 
file and then back to a PLOT3D file will give you an identical file. 


23.4 Using paged meshes and fields 

A program uses paged meshes and fields in nearly the same way as it uses in-core 
meshes and fields. The main difference is how these objects are created. In-core meshes 
and fields can be created by reading a file or by giving a block of data to a constructor. 
Paged meshes and fields should only be created by using functions similar to the ones 
used to read a file for an in-core mesh or field. 
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For example, the FEL.plot3d_read_mesh function will read a PL0T3D file 
into memory and return a FEL_mesh_ptr. The FEL_plot3cLreacLpaged-mesh 
function will open the paged- format file and initialize the paging data structures, and 
then return a FELjmesh.ptr. Both functions take the same arguments. The only 
difference is that the FEL_plot3d_read_pagedmesh function ignores the flags 
argument. See Chapter 22 for more details about these functions. 

Your program can also use file-reading functions that work for both PLOT3D 
files and paged files. Using these generic functions is encouraged since your pro- 
gram will then work with both in-core and paged meshes and fields. These func- 
tions will first look at the flags argument to determine the type of the file. If the 
flags argument specifies a file type (i.e., it is non-zero), that type is used. Other- 
wise, the functions look at the file itself to determine the file type. If the file is a 
PLOT3D file, these functions will read the file and return an in-core mesh or field; 
if it is a paged file, they will open the file and return a paged mesh or field. In gen- 
eral, if the in-core function is FEL_plot3cLreacLxxx, the paged function will be 
named FEL-plot3d_reacLpaged_xxx, and the generic function will be named 
FEL_read_xxx. Table 22.1 shows the correspondences between the generic and 
paged-file functions. 

Another minor difference between in-core and paged meshes and fields is how the 
compute JDOundingJbox and compute_minunax functions are implemented. In- 
core meshes and fields must compute the results at run time. The paged versions of 
these functions instead use values computed when the file was converted to a paged 
file and thus are much faster. Because of this, you should use the FEL functions to 
compute these values instead of writing your own. 


23.5 Controlling memory usage 

If you are concerned about the performance of your program, you will probably want 
to avoid having it use more memory than is installed on the system where it is running. 
You can control how the paging code uses memory by two methods. The first method 
lets you control the total amount of memory used. The second method allows you to 
specify which meshes or fields are more important and should be kept in memory in 
preference to other meshes or fields. 


23,5.1 Pool size 

You can control the amount of memory used by changing the size of the memory pool. 
The paging code uses this pool to hold pages from paged mesh and field files that 
are brought into memory. By default, the maximum pool size is 50% of the system’s 
memory. On most Unix systems, the “system memory size” is the val ue returned from 
the getrlimit system call for the maximu m resident set size. Many shells (sh ^csh, 
etc.) allow this value to be changed with the limit command using the memo ryTa sns 
option. For example, limit memoryuse 100m would set the overall maximum 
memory usage to 100 megabytes and the default paging pool size to 50 megabytes. 
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On Solaris systems, the “system memory size” is determined using a different call that 
returns the actual physical memory size. 

The maximum size of the pool can also be changed with two procedure calls. The 
first, FEL-set_paging_memory .megabytes, sets the size of the pool directly. 
The second call, FEL_set.paging-memory_f raction, changes the fraction of 
the system’s memory that is used; the fraction should range from 0 to 1. Both calls 
take a single float argument. The size of the pool can only be changed before open- 
ing the first paged file. Calls to these two functions after opening a file have no effect. 
The two calls are defined as follows: 

extern void FEL_set_paging_memory_megabytes ( float) ? 

extern void FEL_set_paging_memory_f raction ( float ) ; • 

Memory will only be allocated to the pool as data are read from disk. This means 
that, if the maximum pool size is much larger than the files that are used, the pool will 
grow to at most the size of the files. (The pool can be somewhat larger than the total file 
sizes because partially-filled 8x8x8 cubes of data take up a full page when in memory 
but take only a partial page when on disk.) 

One difference between the paged and standard file reader calls is 
that the same data are loaded in memory only once with the paged 
calls. For example, if you call FEL_read-paged-plot3d-density 

and FEL-read-paged-plot3d_density-rnornentum for the same 
file and use data from both, the density values will only be read once. 
But, if you use the non-paged calls, FEL_read-plot3d_density and 
FEL_read.plot3d-densityjnomentum, the density values will be read 
twice and stored twice. 

23.5.2 Page priority hints 

The second memory control method allows a program to control which file’s pages are 
retained in the memory pool. When the memory pool is full and more data are needed 
from a file, a page that currently contains data is selected and reused. Each mesh or 
field has a priority that determines the length of time before its pages can be reused. 
Increasing or decreasing the priority gives the paging code a hint on which data are 
important or unimportant so that pages containing less important data can be reused 
first. These hints have not yet been used in an application, so they should be regarded 
as experimental. 

For example, priority hints could be used with an unsteady dataset in an interactive 
application that has a concept of the current simulation time. In such an application, 
only the two timesteps that bracket the current simulation time are used. If all the 
timesteps but the current ones are set to low priority, then they will always be the ones 
reused. 

The following priorities can be used: 

• FEL-LOW.P AGE-PRIORITY allows the mesh’s or field’s pages to be deallo- 
cated and reused immediately. 
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• FEL-STANDARD -PAGE-PRIORITY allows the mesh’s or field’s pages to be 
deallocated after they have not been used for a while. This is the default pri- 
ority. 

• FEL -HIGH-PAGE-PRIORITY allows the mesh or field’s pages to be deallocated 
after they have not been used for a longer time than standard-priority pages. 

Priorities are set by using the set member function with one of the following 
keywords: 

• FEL-FIELD-PAGE.PRIORITY changes the priority for both the field (solution) 
data and the mesh data when invoked on a field, and changes the priority of the 
mesh data when invoked on a mesh. 

• FEL-MESH-PAGE.PRIORITY changes the priority of the mesh (if invoked on a 
field, changes priority of the field’s mesh). 

• FEL-SOLUTION-P AGE-PRIORITY changes the priority for the field data when 
invoked on a field, and has no effect when invoked on a mesh. 

The current interface allows you to set the priority for an entire field, an entire 
mesh, or both, or for an individual zone in a mesh. You cannot set the priority for 
individual zones of the data for a multi-zone field. Some examples are below: 

void priori ty_example ( ) 

( 

// code assumes "grid" and "solution* files are 
// multi-zoned 

unsigned int flags = FEL_deduce_mesh_type ( "mesh" ) ; 
FEL_mesh_ptr mesh = FEL_read_mesh ( "mesh" , flags); 
FEL_float_field_ptr field = FEL_read_density (mesh, 

"solution", flags); 

// change priority of mesh for all zones 

mesh->set { FEL_MESH_PAGE_PRIORITY , FEL_LOW_PAGE_PRIORITY) ; 

// change priority of mesh for zone 1 
mesh->get_zone ( 1 ) ->set ( FEL_MESH_PAGE_PRIORITY , 

FEL_LOW_PAGE_PRIORITY) ; 

// change priority of mesh for all zones 

f ield->set { FEL_MESH_PAGE_PRIORITY, FEL_LOW_PAGE_PRIORITY) ; 

// change priority of mesh for zone 1 

f ield->get_mesh ( ) ->get_zone ( 1) ->set (FEL_MESH_PAGE_PRIORITY, 

FEL_LOW_PAGE_PRIORITY) ; 

// change priority of field data for all zones 

f ield->set (FEL_SOLUTION_PAGE_PRIORITY, FEL_LOW_PAGE_PRIORITY) 
// change priority of mesh and field data for all zones 
field->set (FEL_FIELD_PAGE_PRIORITY, FEL_LOW_PAGE_PRIORITY ) ; 


) 



Chapter 24 

The PLOT3D Field Manager 


The PLOT3D field manager is a class hierarchy which can help create and manage 
fields based on PLOT3D data files. Once you create a field manager object, you can 
ask it to create any of over fifty predefined derived fields by name, and it will take care 
of any necessary file reading and derived field construction, returning a field ready for 
use. A function which receives a pointer to a field manager object essentially has been 
sent an entire collection of fields. 


24.1 Constructing an FEL_plot3d_f ield 

The base class of the field manager hierarchy is FEL.plot3d_field. This class 
is an abstract class. Three derived classes create instances of FEL_plot3d.f ield 
for three types of fields: steady fields, time-varying fields, and PLOT3D Q 
fields. Most applications should assign instances of the derived classes to a 
FEL.plot3d_f ield.ptr so that their code is independent of the actual field man- 
ager type. 

24.1.1 Constructing an FEL_steady_plot3d_f ield 

Two constructors are available for constructing an FEL-Steady.plot3d-f ield: 

FEL_steady_plot3d_f ield (char* mesh_file, char* soln_file, 

unsigned flags = 0) ; 

FEL_steady_plot3d_f ield(FEL_mesh_ptr mesh, char* soln_file, 

unsigned flags = 0) ; 

In the first constructor, you supply the names of your PLOT3D mesh file (mesh.f i 1 e) 
and solution file (soln.f ile), and, optionally, flags describing your file formats. In 
the second constructor, you supply a pointer to a preexisting mesh (mesh) instead of 
the name of a mesh file. 
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24.1.2 Constructing a time-varying field manager 

Time-varying field manager objects create time- varying fields in response to requests 
for derived fields. Since manager objects create time-varying fields, the arguments 
to the two constructors for a time-varying field manager are nearly the same as the 
constructors for time- varying fields. Chapter 25 describes time-varying fields and the 
common arguments. 

Three of the field manager constructor’s arguments are different from the time- 
varying field’s constructor. The first argument can be a string specifying the mesh file 
instead of a mesh object. The third argument, the callback function, is different in that 
it does not return a field object. Instead, it fills in the name of the solution file and 
returns a value indicating success or failure. Finally, the field manager has a seventh 
argument giving the file flags for the solution files. 

The definitions for the constructors are: 

const int FEL_plot3d_f ilename_len = 1024; 

typedef int ( *FEL_plot3d_f ilename_callback) 

(int, void*, char [FEL_plot3d_f ilename_len] ) ; 

FEL_f ixed_interval_time_series_plot3d_f ield( 
char* mesh_f ile , 
int num_time__steps , 

FEL_plot3d_f ilename__callback cb, 
void* user_data, 

float physical__timeO , ' • : * 

float physical_timel , 
unsigned flags = 0) ; 

FEL_f ixed_interval_time_series__plot3d_f ield( 

FEL_mesh_ptr mesh, 
int num_time_steps , 

FEL_p 1 o 1 3 d_f i 1 ename_c all back cb, 
void* user_data, 
float physical_timeO , 
float physical_timel , 
unsigned flags = 0) ; 

You can specify time-varying meshes using the second constructor and passing a 
time- varying mesh as the first argument (see Chapter 26). 

An example of a time-varying manager constructor is shown below. The exam- 
ple code reads a steady mesh file called mesh and 20 solution files named solnOO, 
solnOl solnl9. The physical time values for the timesteps start at 0 and in- 

crease by 1 for each timestep. The code also has an example of how to use the user- 
data argument: the first four letters of the filename are passed to the callback function 
using that argument. 
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int soln_callback ( int timestep, void* userdata, 

char filename [FEL_f ixed_interval_time_series_plot3d_f ield 
FEL__plot3d_f ilename_len] ) 

( 

sprintf ( filename , "%s%02d", (char*) userdata, 
timestep) ? 
return 1; 

} 

FEL_plot3d_f ield_ptr make_time_va:r/ing_inanager ( ) 

{ 

unsigned flags = FEL_deduce_mesh_type ( "mesh" ) ; 
FEL_plot3d_f ield_ptr fp = new 

FEL_f ixed__interval_t ime_series_p lot 3d_f ield ( "mesh" , 20 , 
soln_callback, "soln", 0, 1, flags); 
return fp; 


24.1.3 Constructing an FEL_plot3d_q_f ield 

The third type of field manager constructor allows using fields that are not handled 
by the first two cases, such as transformed fields or multizoned fields with mixed 
steady and unsteady zones. This constructor allows an application to create an arbi- 
trary PLOT3D Q field and then to use the field manager to create derived fields from 
it. The constructor simply takes the Q field as its only argument: 

FEL_q_plot3d_f ield(FEL_plot3d_q_f ield_ptr q_field) ; 

This class name can unfortunately be easily confused with another 
type name, FEL_plot3d_q-f ield. The latter type is a synonym for 

FEL_typed_f ield<FEL_plot3d_q>. 


24.2 Creating primitive and derived PLOT3D fields 

Once you have a properly initialized FEL.plot3cLf ield, creating the PLOT3D 
primitive and derived fields is straightforward. Access to the primitive PLOT3D fields 
is provided by five methods: 

FEL_f loat_f ield_ptr get_density_f ield ( ) ; 

FEL_vector3f_f ield_ptr get_momentum_f ield( ) ; 

FEL_f loat_field_ptr get_energy_f ield{ ) ; 

FEL_p 1 o t 3 d_dens i ty_momen t um_f i e 1 d_p t r 
get_density_momentum_f ield{ ) ; 

FEL_plot3d_q_f ield_ptr get_q_f ield ( ) ; 
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The first three methods return fields representing the momentum variable and two ther- 
modynamic state variables which make up the PLOT3D solution output. The latter 
two methods return “conglomerate” fields, the last one containing the entire PLOT3D 
solution vector. 

PLOT3D derived fields are obtained via two functions: 

FEL_f loat_f ield_ptr make_f loat_f ield (FEL_f loat_f ield_enum) ; 

FEL_vector3f_f ield_ptr 

make_vector3 f_f ield { FEL__vec tor 3 f_f ield_enum) ; 

These functions take a single argument indicating the desired field and they return a 
pointer to a field of the appropriate type. The enums representing the supported fields 
are shown in Table 24.1. The identity of the fields produced by the various enums 
should be evident from their names; precise definitions can be found in [WBPE92] or 
in the PEL code itself (in the file FEL.pl ot3cLmap-f unctions . C). The enums 
correspond exactly to the original PLOT3D “function numbers” ([WBPE92]), so you 
can use those directly if you prefer. If you need a derived field which isn’t predefined, 
you can always create one using some combination of the primitive and predefined 
fields — see Chapter 19. 


FEL.density = 100 

FEL-stagnat ion-density * 102 

FEL.pressure * 110 

FEL-st agnation_pressure * 112 

FEL.pressure.coeff icient - 114 

FEL.pitot-pressure = 116 

FEL.dynamic .pressure = 118 

FEL-normali zed-temperature = 121 

FEL_normali zed_stagnat ion-temperature - 123 

FEL-normali zed-enthalpy - 131 

FEL-normali zed-stagnation enthalpy * 133 

FEL-normali zed-internal energy = 141 

FEL-normali zed-stagnation energy = 143 

FEL-normali zed-kinetic.energy = 145 

FEL.v.velocity - 151 

FEL-velocity-magni tude = 153 

FEL.speed.of .sound = 155 

FEL-divergence.of.velocity * 158 

FEL.y .momentum = 161 

FEL.energy =163 

FEL.entropy.sl = 171 

FEL.y .component. of -vorticity = 181 

FEL.vorticity-magnitude = 183 

FEL.veloci ty.crcss.vorticityunagni tude = 185 

FEL.pressure.gradientunagni tude = 192 

FEL-shock = 400 


FEL-normali zed-density = 101 

FELjnormalizecLstagnation-density = 103 

FEL-normali zed-pressure = 111 

FEL-normalized-stagnation-pressure = 113 

FEL-stagnation-pressure.coef f icient * 115 

FEL-pi tot-pressure-ratio = 117 

FEL-temperature =120 

FEL_st agnation- temper a ture = 122 

FEL-enthalpy = 130 

FEL-stagnat ion-enthalpy = 132 

FEL-intema 1-energy = 140 

FEL-stagnat ion-energy = 14 2 

FEL-kinetic.energy = 144 

FEL_u-velocity = 150 

FEL-W-velpcity = 152 

FELjnach-n umber = 154 

FEL-cross flow velocity = 156 

FELjc-momentum = 160 

FEL.z unomen turn = 162 

FEL.entropy = 170 

FEL-X-component-of.vorticity = 180 

FEL-Z-Component.of.vorticity = 132 

FEL-Swirl = 184 

FEL_helicity = 186 

FEL-density.gradient magnitude = 193 
FEL-f iltered-shock = 401 


FEL.velocity = 200 
FEL unomen turn = 202 
FEL-velocity-cross.vorticity = 204 
FEL_density_gradient = 211 


FEL-vorticity = 201 
FEL_perturbation-velocity = 203 
FEL.pressur e-gradient = 210 


Table 24. 1 : Scalar (above the line) and vector (below the line) derived fields predefined 
in PLOT3D [WBPE92] and supported by FEL.plot3d.f ield. The PEL enums are 
given, along with the PLOT3D “function number”. 
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AH" the derived fields created by the FEL.plot3d-f ield are of type 
FEL_map-then_interpolate_derived_f ield* ( ) . Any required differential 
operator fields are second-order if the FEL_plot3d-f i eld’s mesh is structured, and 
first-order if the mesh is unstructured. Future versions of the FEL.plot3d-f ield 
manager may permit the client to specify the types of derived fields and differential 
operator fields to be used in the construction of an individual field. 


24.3 How the field manager works 

The PLOT3D field manager works by reading core or paged fields and return- 
ing derived fields from these underlying fields. When it creates a core or paged 
field, it saves a pointer to it and when possible reuses the field for later de- 
rived fields. For example, if you call make.f loat.f ield (FEL_pressure) and 
make.f loat.f ield ( FEL.temperature ) , they will share the underlying Q field. 

This sharing of multiple derived fields only happens when the the same field man- 
ager object is used to create all of the derived fields. An application that creates a field 
manager object, creates a derived field, destroys the field manager, and then creates a 
new field manager and a second derived field will read the solution twice. 

Memory usage with the field manager is not always optimal when the field manager 
is using PLOT3D solution files. When a field manager is given a PLOT3D solution file, 
in most cases it reads the entire file and creates a Q field when the first derived field is 
requested. This is not efficient if only part of the solution field is used, which would 
happen if make.f loat-f ield (FEL.density) is called. However, reading all of 
the data allows multiple derived fields to share the solution data when the derived fields 
need different but overlapping portions of the solution file. 

The primitive PLOT3D field functions (the get.* -field functions) for steady 
and time-varying field managers operate differently. Each reads a portion of the solu- 
tion the first time it is called. If all file primitive field functions are called, the density 
and energy solution data will be read twice, and the momentum data will be read three 
times. 

If you first create a primitive PLOT3D field and then create a derived 
field, the field manager may use the primitive field instead of first creat- 
ing a Q field. For example, if get-density-momentum-field and then 
make_vector3f.f ield ( FEL.velocity) are called, the velocity field will be 
derived from the density-momentum field created in the first call. If the density- 
momentum field had not existed, the velocity field would be derived from a Q field. 
This behavior depends on which primitive PLOT3D fields are needed for each derived 
field. See the source file FEL-plot3d_f ield . C for details. 

Memory usage with paged files is not an issue. The paging code insures that the 
data needed is read and stored only once, even if it is used for multiple primitive 
PLOT3D and derived fields. Memory usage is also not an issue when the field manager 
is given a PLOT3D Q field. In this case, the primitive and derived field functions return 
fields derived from the Q field handed to the field manager constructor. 
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24.4 Miscellaneous FEL_plot3d_f ield methods 

FEL.plot 3d_f ield provides a get_mesh ( ) and a generic set ( ) method as a 
convenience. This code: 

FEL_plot3d_field_ptr plot3d_manager; 

FEL_mesh_ptr mesh = plot3d_manager->get_mesh ( ) ; 

returns the mesh associated with plot3djnanager. This mesh is associated with 
any of the fields created by plot3d_manager. 

FEL_plot3d_f ield: : 

set (const FEL_set_keyword_enum keyword, int value); 

passes the requested options to the mesh and primitive fields of the 
FEL-plot3d_f ield object. Note that since all fields created by the 

FEL_plot3d.f ield share the same mesh and primitive fields, the set ( ) 
call on a FEL.plot3d.f ield object affects all fields which have been created, 
and even fields which will be created, by the FEL_plot3d_f ield. Furthermore, 
set ( ) calls on any field created by the FEL-plot3d_f ield will have the same 
wide-reaching effect. At present, this interdependence of a group of fields sharing a 
common mesh is a limitation of FEL which will be resolved in a future release. In the 
meantime, be aware of potentially unwanted side effects of the set ( ) calls. 


24.5 An example 

# include "FEL.h" 

int main (int argc, char** argv) 

{ 

// mesh file, solution file in argvfl] and argv[2] 

FEL_plot3d_f ield_ptr plot3d_manager = - 

new FEL_steady_plot3d_f ield (argvfl] , argv [2] ) ; 

FEL_f loat_f ield_ptr helicity_f ield = 
plot3d_manager->make_f loat_f ield (FEL_helicity ) ; 

FEL_vector3f_f ield_ptr grad_pressure_field = 
plot3d_manager->make_vector3f_field(FEL_pressure_gradient) 

// gradjpressure^f ield is ready for at_calls, or whatever 

return 0; - = 

} 
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This example first creates a FEL.steady.plot3d_f ield using the mesh and so- 
lution files specified in the first two arguments. It then creates a float field containing 
helicity and a vector field containing the gradient of pressure. 


24.6 PLOT3D derived field “convenience functions” 

The FEL_plot3d_f ield takes care of file reading and caching of primitive fields by 
maintaining internal state, but in response to an incoming FEL_f loat.f ield_enum 
or FEL-vector3f.f ield_enum the manager creates the requested derived field by 
calling one of a number of external "convenience functions”. You can call these func- 
tions directly, whether or not you’ve created an FEL_plot3d_f ield. All of the 
functions take a pointer to an FEL-plot3d_q_f ield as their only required argument 
and, as usual, an optional character string name. A few of the functions have over- 
loaded versions which accept pointers to fields comprising a subset of the PLOT3D 
solution vector. The functions return a pointer to either an FEL-f loat.f ield or an 
FEL.vector3f -field. For example: 

FEL_plot3d_q_f ield_ptr q_field; 

FEL_f loat_f ield_ptr sp_field = 

FEL_plot3d_make_stagnation_pressure_f ield (q_f ield) ; 

Although this looks very similar to a derived field constructor invocation, this wrapper 
function calls any required constructors on your behalf, so you don’t want a new after 
the 

Tables (24.2 and 24.3) present the available PLOT3D derived field convenience 
functions. These functions correspond one-to-one with the FEL enums shown in 
Table 24.1, modulo an obvious naming convention. The first column in the table 
lists the function. The second column lists the arguments accepted by the func- 
tions. Many functions are overloaded so multiple arguments types are accepted. 
Table 24.4 lists the abbreviations for the argument types. See the header file 
FEL-plot3d_derived_f ield.h for more details. 

The default name of the created field is the name of the function called with the 
FEL_ and make, deleted. For example, FEL.plot3d-make-density-f ield cre- 
ates a field with a default name of FEL _plot3d.density -field You can override the de- 
fault name by specifying the name as the second argument (the third argument for the 
FEL-plot3d-make.velocity_f ield function that accepts separate density and 
momentum fields). 


24.7 PLOT3D derived field inventory arrays 

An application may want to present all PLOT3D derived fields to a user as a selectable 
list. Since the fields are predefined and lazily evaluated, they can all be created ahead 
of time with minimal time and storage overhead and, when selected, they can be easily 
conjured forth without any need for a complicated input routine and interpreter or the 
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Function Name 


FEL.pl ot3d_make_density-f ield 


FEL.pl ot3d.make-normali zed-density-field 


FEL-plot 3d-make_stagnati on-dens lty-field 


FEL-plot3djnake-normalized-stagnation..density-f ield 


FEL-plot3d-make-pressure.f leld 


FEL.pl o 1 3d-make_normal i z ed-pressure.f i eld 


FEL-plot3d-make_stagnation.pressure-field 


FEL-plot 3 d_make-normalized_stagnation.pressure-field 


FEL.pl ot3d.niake.pressure_coe£ f icient-f leld 


FEL-p 1 o 1 3 d-make-s t agna t i on_pressure-coef f icient-f leld 


FEL_plot3djnake_pi tot pressure-field 


FELplot3d-make -pi tot-pressure -ratio-field 


FEL.pl ot3cLmake -dynamic-pressure -field 


FELplot3d_make-temperature.f leld 


FEL-plot3d-makeJiormali zed-temperature-field 


FEL_plot3d_make-stagnat ion-temper aturepield 


Arguments 


Q 


Q,D 


Q 


Q 


Q 


Q 


Q 


Q 


Q 


Q 


Q 


Q 


Q, DM 



FEL-plot 3d-make.enthalpy-f ield 

Q 

FEL-plot3d-make_normali zed.enthalpy-field 

Q 

FEL.plot3d_make.st agnation-enthalpy-field 

Q 


FEL-p 1 ot 3d-make-normali zed-stagnat ion-enthalpy -field 


FEL.p lot 3d _make-internal .energy-field 


FEL.p 1 o 1 3 d _ma ke -no rma 1 i z e d-i n t ema 1 -energy-field 


FEL.p lot3d-make-$t agnation-energy-field 


FEL_plot3d-make.normali zed-5 tagnation.energy_f leld 


FEL.p lot 3d jnakeJcinetic.energy.f ield 


FEL.p 1 o 1 3 djnake -normal i zed-kinet ic.energy_field 


FEL_plot3djnake_u.velocity-f ield 


FEL-plot3d-make-V-velocity.fi eld 


FEL-plot 3d_make-W.velocity.f ield 


FEL.plot 3d_make.veloci ty_magni tude.f ield 


FEL_plot3d_make-mach-n umber -field 


FEL-plot 3 d-make-speed.o f -sound.f i eld 


FEL-plot3d-make.cross-f low -velocity-field 


FEL-plot3d-make.divergence.of .velocity-field 


FEL.plot 3d-make_x -momentum-field 


FEL.plot 3d-make-y-momentum_f ield 


FEL-plot 3 d-make.z .momentum-field 


FEL.plot 3d_make.energy-f ield 


FEL.plot 3d-make.entropy-f ield 


FEL-plot 3d_make-entropy_sl-f ield 


FEL-plot 3dmake-x-component-o£-vorticity-f ield 


FEL-plot 3d-make-y-Component-0 f -vortici ty.f ield 


FEL.plot 3d-m.ake_z -component -O f-vorticity-f ield 


FEL-plot 3d-make-vorticity-magnitude-f ield 


FEL-plot 3d-make-swirl.f ield 

FEL-plot 3d-make-velocity.cros5-vorticity-magnitude.f ield 
FEL_plot 3d_make_helicity-f ield 


FEL-plot 3 d_make-pressure-gradient_magnitude-f ield 


FEL_plot 3d_make_dens i ty.gr adi ent-magni tude.f ield 



Q 


Q 


Q, DM 


Q, DM 


Q, M 


Q, M 


Q, M 


Q 


Q 


Q 


Q, DM 


Q, DM 


Q, DM 


Q, DM 


Q, DM 
Q, DM 
Q, DM, V 


Q 


Q.D 


Table 24,2: Derived field convenience functions that return float fields. 
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Function Name 

Arguments 

FEL.ploc3d_make-velocity_f ield 

Q, DM, D+M 

FEL.pl o 1 3 djnake.vor t i c i ty.f i e 1 d 

Q, DM 

F E L.p 1 o 1 3 cLmak e jnomen t uin_ f i e 1 d 

Q 

FEL.plot3djnake^)erturbation_veloci ty_f ield 

Q, DM 

FEL-plot3djnake-velocity-cross-vorticity.f ield 

Q, DM 

FEL_plot3djnake-pressure_gradient-f ield 

Q 

FEL-plot3d-n\ake_density_gradient_f ield 

Q.D 


Table 24.3: Derived field convenience functions that return vector fields. 


Abbreviation 

Argument Description 

D 

Density field 

DM 

Density -momentum field 

D+M 

Two separate arguments: a density field and a 
momentum field 

M 

Momentum field 

Q 

Q field 

V 

Velocity field 


Table 24.4: Abbreviation key for the derived field function tables. 


like. As an aid for this kind of usage, FEL provides two arrays containing all the FEL 
PLOT3D enums — one with all the float field enums and one with all the vector field 
enums. A sentinel (“0”) marks the end of each array. 

The “inventory arrays” are found in FEL_plot3d_f ield.h: 

FEL_f loat_f ield_enum FEL_plot3d_f loat_f ields [ ] ; 

FEL_vector3f_f ield_enum FEL_plot3d_vector3f_f ields [] ; 

By traversing these arrays, one can sequentially construct all PLOT3D derived 
fields and, by calling get_name ( ) on each field, one can accumulate a list of canon- 
ical field names to be presented to the user. If you are listing the names, the field 
manager should be created using a FEL_plot3d_q_f ield instead of actual files. 

If you use actual files, the names may vary from the canonical plot 3 d_xxc.fi eld 
form. Here is a short demonstration showing how to list all the derived scalar fields 
using FEL_plot3d_f loat.f ields [ ] : 

# include " FEL . h " 

char** get_derived_f loat_names ( int* nsf ields) 

{ 

// make dummy field manager using dummy mesh and field 
F EL_me s h__p t r mesh = new FEL_regular_mesh (4 , 4, 4) ; 
FEL_plot3d_q_f ield_ptr q_field = 

new FEL_constant_f ield<FEL_plot3d_q> (mesh, FEL_j?lot3d_q ( ) ) 
FEL_solution_globals sg; 

sg . set_alpha ( 0 ) ; sg . set_f ree_stream_mach ( 0 ) ; 
sg . set_t ime ( 0 ) ; sg . set_reynolds_number ( 0 ) ; 
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q_f ield->set_solution_globals (sg ) ; 

FEL_plot3d_f ield_ptr plot3d_field = 
new FEL_q__plot3d_field(ci_field); 
assert (plot3d_f ieldl =NULL) ; 

*nsf ields = 0; 

// no sizeof on unspecified dim array 

while ( FEL_plot3d__float_f ields [ ( *ns fields) ++] != 0) ; 

-- ( *nsf ields) ; 

FEL_f loat_f ield__ptr* scalar_f ields = 
new FEL_float_field_ptr [*nsf ields] ; 
char** scalar_names = new char* [ *nsf ields] ; 
for (int i=0; i<*nsfields; ++i) { 

FEL_f 1 oa t_f i e 1 d_p t r f = 

plot3d_f ield->make_f loat_f ield 
(FEL_plot3d_float_f ields [i] ) ; 
assert ( f ! =NULL) ; 
const char* n = f ->get_name ( ) ; 

scalar_names [ i ] = strcpy(new char [ strlen (n) +1] , n) ; 
f =NULL»; 

} 

// make names without "plot3d tt # " field” , and 
char t [256 ] ; 
char* s; 

for ( i = 0 ; i< *ns fields ; ++i ) { 
t[0]='\0 ' ; 

for ( s=scalar_names [ i 3 ; ( s=strtok (s , ) ! =NULL; s=NTJLL) 

if ( strcmp ( s , "plot3d H ) && strcmp ( s , " field" ) ) { 

s treat ( t, s) ; 
s treat ( t , " " ) ; 

} 

t [strlen(t) -1] = '\0'; 

// t can't be longer than original 
strepy { scalar_names [ i ] , t) ; 

} 

return scalar_names ; 


int main() 

{ 

- int nsf ields; 
char** scalar_names = 

get_derived_f loat_names (Sens fields) ; 
cout << nsfields « 

n scalar fields ready and waiting in scalar_f ields [] " 
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<< endl; 

for (int i=0;i<nsfields;++i) 

cout << scalar_names [i] « endl; 
return 0; 


} 
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Chapter 25 

Time-Varying Fields 


Fields that vary with time, also known as unsteady fields, are often represented by a 
series of time steps, where each time step represents the field at some instant in time. 
The time steps may or may not be at fixed intervals during a simulation, though in 
many cases a fixed interval is used. FEL represents unsteady fields via the classes 
FEL_time.series-f ield and FEL-f ixed.interval-time-series-f ield. 
The interface for time-series fields inherits from the standard field interface defined in 
FEL-f ield and FEL-typed.f ield. Just as with steady fields, one can make “at” 
calls to request data using an FEL_phys-pos, FEL-cell, FEL-vertex.cell or 
FEL.structured.pos argument. All the argument types contain an FEL-time 
data member. For steady data, the time is ignored; for unsteady data, it is essential 
that time be specified. Given an argument containing a time value, time-series fields 
produce field values by accessing data from the appropriate time step. If necessary, 
time-series fields can also do temporal interpolation (described later in this chapter) to 
produce values at a time intermediate to the time steps. 

In general, a time-varying field with node type T can be used anywhere a field of 
node type T can be used; for instance, one can build derived fields on top of an un- 
steady field. There are a few cases where, due to precomputation or caching, unsteady 
fields are not interchangeable with steady fields. Recall that the member function 
get .eager.f ield provides a way for the user to get a new field where every node is 
eagerly evaluated. With a time-varying field, get.eager.f ield essentially returns 
a snapshot in time, i.e., a steady field. Thus if one wants to call get.eager.f ield, 
then the final argument, specifying a time value, is no longer optional. It is a run-time 
error if one calls get_eager_f ield on an unsteady field without providing a time. 

A second case where unsteady fields require special treatment occurs when 
working with the FEL caching fields. The fields FEL-cachecLf ield and 
FEL Jnash.cachecLf ield provide an option for users who want the computational 
frugality of caching results, without having to pay the costs of eagerly precomputing 
values over a whole field. Unfortunately, the FEL caching fields are not designed to 
work with time-varying fields; caching fields ignore time. Thus, with caching fields 
as they are currently implemented, a cached time-varying field would result in better 
performance, but in general wrong answers. To prevent what could be anr insidious 
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problem, the cached field constructors check the type of field they are caching (at run- 
time) and fail if the argument is not steady. 

FEL currently supports time-series data where there is a fixed interval in time be- 
tween each time step. The library is designed so that adding support for data that are 
not temporally spaced at fixed intervals should not be too difficult in the future. The 
library supports both time- varying fields and time-varying meshes. The time-varying 
mesh support is described in the next chapter. The topology of the mesh is currently 
required to be constant over time, in other words, meshes (and the fields based on those 
meshes) that adapt over time are not cunently supported. 

An individual time step in an unsteady data set may be very large, the set of all 
the time steps may be hundreds of times larger. For many data sets it is not feasible 
to load every time step into memory. FEL time-series fields support the automated 
management of a subset of the time steps in memory using a working set. We describe 
the FEL-time.series-f ield working set control interface next. 


25.1 Working sets and callbacks 

The FEL_time_series,f ield class manages a working set of field objects, each 
object corresponding to a time step. When the user requests data at some point within 
the field represented by the FEL.t ime.series,f ield instance, the time-series field 
must check to see if the necessary time steps are currently in the working set. If so, 
then the time-series field requests data from the time steps, does temporal interpolation 
if necessary, and produces the result. If a needed time step is not in the working set, 
and the working set is not full, then the time step is acquired using the user-provided 
callback function. The callback function is described below. If the working set is full, 
then the field replaces a time step using a least recently used policy to choose the field 
to be replaced. 

For most users, the working set management is completely automatic, in other 
words there is no need for the user to manually load and unload data. For those who do 
want to directly control which time steps are loaded, the FEL-time_series_f ield 
class provides an interface where one can dictate the working set size and contents. 
The working set management member functions are: 

void set_working_set_size ( int n) ; 
bool loaddnt time_step} ; 
bool load_all(); 
void unload(int time_step) ; 
void unload__all ( ) ; 

void unload_least_recently_used( ) ; 

void set_verbose (bool ) ; 

void show__working_set (ostream& s) ; 

The function set.working-set.size allows the user to control the size of the 
working set used by the field. The default working set size is 5. The function load 
allows the user to make the field load a particular time step. The function load_all 
resizes the working set to the total number of time steps and loads each one. The 
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loacLall function is handy when an application can afford to load every time step. 
The unload functions are relatively self-explanatory. The set_verbose function 
can be used to make the field output updates when time steps are loaded and unloaded. 
Watching the output from verbose mode can be educational for those not familiar with 
the use of working sets. 

Users working with applications that are multi-threaded need to be careful with 
time-series fields as they are currently implemented. The routines that manage the 
working set do not have critical section protection. This means that, if multiple threads 
make accesses that cause the working set to change, then there is a chance that the 
working set data structure will become inconsistent. To avoid this problem the user 
should ensure that the working set does not change while in multi-threaded code. The 
simplest way to ensure no change is to load all the time steps initially. The load-all 
call is handy for this purpose. If loading all the time steps is not an option, then the user 
must make sure that working set changes occur only in single-thread mode. Providing 
critical section protection for the working set is on the list of future enhancements to 
FEL. 

FEL-tin\e_series_f ield relies on a callback function provided by the user at 
field construction time to produce the field corresponding to a particular time step when 
needed. For example, a callback returning a density field might look like: 


unsigned flags = FEL_PLOT3D_3D_WHOLE; 

const char* f ile__names [ ] = ("filel", H file2", "file3", "file4 

FEL_f loat_f ield_ptr my_callback (FEL_mesh_pt r m, int time_step 

void* ) 

{ 

return FEL_read_density (m, file_names [timers tep] , flags); 

} 

The callback function takes three arguments: the mesh rn that the return field should 
be based on, a time step, and a void* pointer to “client data”, i.e., data provided by 
the user when the time-series field is constructed that gets handed back to the callback 
function. The client data, for instance, could be a pointer to a structure that contains 
parameters such as the file names and file reader flags, parameters that are globals in 
the example above. The callback should return the field for the given time step. The 
callback can return NULL to indicate some type of failure, if a file could not be read, 
for example. 

Callbacks give the user a great deal of flexibility in complying with requests for 
time step data. For example, the callback could construct and return a derived field. 
Another possibility would be for the callback to construct a new FEL.core.f ield 
using a buffer already in memory, such as the buffer used by a simulation. Yet another 
possibility would be for the callback to construct some type of analytic field that could 
be used for testing. 
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25.2 Time representations and conversions 

The objects in FEL, such as FEL_phys_pos, that contain a time value use the type 
FEL-time. FEL_time, as described in Chapter 4, can represent three types of time: 
physical, computational or integer computational (time step). Usually it is not neces- 
sary to convert from one representation to another, since the user can choose to work 
in any of the three representations. If the user does need to convert, then the library 
provides the function: 

int convert_time (const FEL_time& from, 

FEL_time_representation_enum to_representation, 
FEL_time* to) const; 

The routine converts from from to to, given the desired to_representation. 


25.3 Temporal interpolation 

Given a time value that does fall on a time step, time-series fields do temporal inter- 
polation. FEL currently supports linear interpolation between the pair of time steps 
bracketing a particular time. No temporal interpolation will take place if the user is 
working in the integer time step representation. Temporal interpolation is also sup- 
pressed if working in physical or computational time, and a given time value implies 
no fractional part between time steps. 

In the case where both temporal and spatial interpolation are necessary, for instance, 
to produce a field value at a physical position in an unsteady field, temporal interpola- 
tion occurs first. For instance, consider a hexahedral structured mesh with simplicial 
decomposition turned off. Given an at _phys_pos call with a point p, FEL will locate 
the hexahedron containing p, temporally interpolate to get the field values at the time 
specified in p and then spatially interpolate to get the final result. 

25.4 A time- varying field example 

To illustrate the construction and use of a time-varying field, we present a small ex- 
ample. The callback and globals used for this example are the same as used for the 
example above (my .callback). See also the example program included as part of 
the FEL primer: primer.l3a . C. 

#include "FEL.h" 

int n_time_steps = 4; 

float physical_time_0 = 12000.0; 

float physical_f ixed_interval = 1.0; 

int main ( ) 

{ 

int res; 

flags = FEL_deduce_mesh_type(argv[l] ) ; 


25 A. A TIME- VARYING FIELD EXAMPLE 


153 


FEL_mesh_ptr mesh = FEL_read_mesh (argv[l] , flags); 
assert (mesh != NULL) ; 

FEL_f loat_f ield_ptr field = 

new FEL_f ixed_interval_time_series_f loat_f ield ( 
mesh, 

n_time_steps # 
my^callback, 

NULL, 

physical_time_0 , 
physical_f ixed_interval) ; 

// find an arbitrary physical point inside mesh 
FEL_cell cell; 

FEL_phys_pos phys_pos; 

res = mesh->int_to_cell (mesh->card(3 ) / 2, 

3, &cell ) ; 

assert (res == 1) ; 

res = mesh->centroid__of_cell (cell , &phys_pos) ; 
assert (res == 1 ) ; 

phys_pos . set_physical_time (physical_time_0) ; 
float f; 

int res = f ield->at_phys_pos (phys_pos , &f ) ; 
assert (res == 1) ; 

} 

The mapping between physical time and time steps is specified by the physi- 
cal_t ime_0 and physical.f ixed.interval arguments to the field constructor. 
For a time step t, the corresponding physical time p would be: 

p = physical_time_0 -I- 1 * physical-fixed-interval. 

When the user works in physical time, FEL_time_series-f ield solves for 
(floating-point) computational time t using the same equation. If t has no fractional 
part, then temporal interpolation can be avoided. If t does have a fractional part, then 
the bracketing time steps are equal to the floor and ceiling of t. 


154 


CHAPTER 25. TIME-VARYING FIELDS 



Chapter 26 


Time- Varying Meshes 


It is not difficult to imagine unsteady data sets where not only the field values change 
with time but also the mesh geometry. For instance, the simulation of a helicopter or 
windmill typically involves rotating blades, and often there are meshes that move with 
the blades. FEL currently supports time-varying meshes using the time-varying field 
mechanism described in the previous chapter. The time-varying support is currently 
limited to curvilinear meshes. Once the time-varying mesh is constructed, it can be 
used just as any other mesh. As with time- varying fields, the biggest difference when 
using an unsteady mesh is that the arguments to mesh calls (e.g., the FEL.cell passed, 
to coord.inates_at.cel 1) must have their time value set. 

Unfortunately, to construct a time-varying mesh, one must follow a somewhat cir- 
cuitous route. A curvilinear mesh can be made to be unsteady by constructing the 
mesh with a time-series field. The field has node type FEL_vector3f _and.int, 
where the vector part represents coordinates, and the int part represents an IBLANK 
value. A curvilinear mesh constructed with a FEL.vector3f .ancLint field con- 
sults the field whenever coordinates or IBLANK data are needed. If the field provided 
at mesh construction time is a time-varying field, then presto, one has a time- varying 
mesh. 

In the following sections we walk through the construction of an unsteady (single- 
zone) curvilinear mesh and discuss what would need to be done in the multi-zone case. 
Building time- varying meshes is currently one of the more challenging things that one 
can do in FEL, and the reader is warned that more than one pass over this chapter may 
be necessary in order for everything to make sense. In the future our plan is to factor 
out the time-series mechanism from the time-series field class, so that it can be used for 
both meshes and fields, without some of the gymnastics described below. In the mean 
time, w'e hope that the following example is illuminating. 

In the example we omit the construction of the solution field. Typically a data set 
with an unsteady mesh also has unsteady solution data. See the previous chapter for the 
details on building time-varying fields. See also the program primer JL3b.C from 
the FEL primer for an example where a time- varying mesh is constructed. 
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26,1 Single-zone time-varying meshes 

To construct an unsteady curvilinear mesh, follow these steps: 

1. Assume a few globals, e.g.: 


const int n_time_steps = 5; 
const char* mesh_names [n_time_steps] = 
{"mO", "ml", "m2", "m3", "m4"> ; 
unsigned f ile_reader_flags,- 
const float physical_time_0 * 0.0; 
const float physical_f ixed_interval = 1.0; 


2. Provide a callback to be used by the FEL.vector3f_andJ.nt 
use an adapter to make a mesh behave like a field, e.g.: 


field, where we 






— — * * *- — *■ ■‘■'-■‘■u i 

my _mesh_ callback (FEL_mesh_ptr, int time_step, void*} { 

F EL_me s h_p t r mesh; 

mesh = FEL_ read_mesh(mesh_names[time_stepj , file_reader flags 

if (mesh == NULL) return NULL; ~ 5 

return new FEL_mesh_as_vector3f_and_int_f ield (mesh) ; 


3. Determine the structured dimensions of the mesh, e.g.: 
int res; 

FEL_mesh_ptr mesh; 

char* mesh_name = "my^mesh"; 

= FEL - d educe_mesh_type(mesh_name) ; 
assert (file_reader_f lags != 0); 

FEL_vector3i dim; 

res - TEL-. read— mesh_zone_ dimensions (mesh_name, 

f ile_reader_ f lags , 0, &dim) 

4. Build a structured mesh: 


mesh - new FEL_structured_mesh(dim[ 0 ] , dim(l] , dim[ 2 ]) ; 

5. Build a time-series field for the coordinates and IBLANK, e.g.: 

FEL_vector3f_and— int-f ield_ptr xyzi_field = 
new FEL_fixed_interval_time_series_vector 3 f_and int field( 

mesh, “ ~~ ' 

n_time_steps, 

my_me sh_ca 1 1 ba ck , 

NULL, 

physical_time_ 0 , 
physical_f ixed_interval) ; 

6. Build the unsteady curvilinear mesh, e.g.: 
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mesh = 

new FEL_curvilinear_mesh_xyzi_f ield_layout (dim[0] , dim[l] , 

dim [ 2 ] , xyzi_field) 

26.2 Multi-zone time- varying meshes 

To construct an FEL _mul t i jnesh, one should follow the same pattern as above, once 
for each zone. If the user callback reads data from a multi-zone mesh file, then typically 
the callback should just read data for a particular zone. If only some zones vary with 
time then a time-series mechanism need only be built for those zones. Thus, if an appli- 
cation has some information about which zones are steady and which are not, there is 
an opportunity for optimization. By not building time- varying objects for steady zones, 
an application should get better performance, since unnecessary temporal interpolation 
and time step loading will be avoided. There should also be memory savings, since the 
mesh for a particular zone will only occur once in memory. For meshes that contain 
PLOT3D IBLANK information, one should be sure that if a zone is assumed to be 
steady, then the IBLANK and coordinate data both should not vary with time. 
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Appendix A 

Glossary 


abstract class An abstract class defines interface and implementation which are in- 
herited by derived classes. Abstract classes in C++ cannot be instantiated. See 
also concrete class. 

abstract factory See factory. 

adapter An adapter provides an alternate interface for a class. For example, in FEL 
the class FEL-mesh.as-f ield<T> serves as an adapter that provides a field 
interface for a mesh object. 

adjacent cells Given n-complex C, by convention there exists one (null) (n + l)-cell 
of which every n-cel! of C is a face; likewise there exists one (null) (n — l)-cell 
which is the face of every vertex. Distinct Ar-cells c and d (for 0 < k < n) are 
then said to be adjacent if: 

(i) there exists some ( k — l)-cell of C that is a face of both c and cf, and 

(ii) there exists some ( k + l)-cell of which each of c and d is a face. 

For example, two hexahedra are adjacent if they share a quadrilateral face. Two 
vertices are adjacent if they are the endpoints of a common edge. 

affine combination Let P = {pi, ...,P*} be a finite set of points in R d . A point x is 
an affine combination of P if: 

k 

x = 

i=l 

where A* = 1. See also linear combination. 

affinely independent A finite set P oik {joints is affinely independent if there is not a 
point pi e P such that pi is an affine combination of P - {pj}. 

block A block is an alternate name for zone. 
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Cartesian grid A Cartesian grid is a grid where the cells are aligned with the coor- 
dinate axes. The grid looks similar to the output from a quad-tree or oct-tree 
decomposition, though Cartesian grids are not necessarily oct-trees. The term 
"Cartesian grid" is also used by some as a synonym for uniform grid. 

cell-centered field A cell-centered field is a field where there is a node associated with 
each cell in the mesh. Typically the cells are hexahedraor tetrahedra. 

CFD Computational Fluid Dynamics. 

cell a k-cell of a topological space T is a subspace of T whose interior is homeomor- 
phic to R* and whose boundary is non-null. Less formally, typically when one 
speaks of cells in CFD one means hexahedra or tetrahedra (i.e., 3-dimensional 
cells or 3-cells), though there are also 0-cells (vertices), 1-cells (edges) and 2- 
celis (polygons). 

cell complex A cell complex of a topological space T is a finite collection C of cells 
of T such that: 

(i) the relative interiors of cells of C are pairwise disjoint, 

(ii) for each cell c € C, the boundary of cell c is the union of elements of <7, 

(iii) if c, d € C and c H d £ 0, then c fi d is the union of elements of C. 

CGNS CONS (Complex Geometry Navier Stokes) is a file format currently under 
development for storing CFD data. 

Chimera Chimera, a composite monster from Greek mythology, is used in CFD when 
speaking of multi-zone grids. 

Chimera file A Chimera file is an auxiliary file to the main multi-zone grid description 
file. The Chimera file contains definitive information on how the solver should 
handle regions where zones intersect. 

class A class is a basic concept in object-oriented programming. A class encapsulates 
both data and member functions for operating on the data. 

closure The closure of a cell c in a cell complex C consists of all the faces of c. 

computational space A computational space is a space defined in terms of a partic- 
ular discretization used to decompose a physical space. One typical computa- 
tional space corresponds to a hexahedral curvilinear (structured) mesh in physi- 
cal space: in computational space each physical space hexahedron corresponds 
to a unit cube. Each vertex in such a computational space can be indexed as if in 
a 3-dimensional array, and the indices are typically labeled i, j, and k. 

concrete class A concrete class is a class that can be instantiated. See also abstract 
class. 

curvilinear mesh A curvilinear mesh is the most general form of structured mesh. 
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demand-driven evaluation Demand-driven evaluation is a technique where compu- 
tations are done only when requested, rather than in advance. 

derived field A derived field is a field defined in terms of one or more other fields. 
For example, a velocity derived field can be defined in terms of momentum and 
density fields. 

dynamic binding The run-time association of a request to an object and one of its 
operations. In C++, virtual functions are dynamically bound. 

dynamic casting Dynamic casting is a new feature in the C++ standard designed to 
support safe casts down or across an inheritance hierarchy. With dynamic casts, 
the success of the cast of a pointer can be verified at run-time by testing whether 
the result is non-NULL. For example, given a class B and a class D derived from 
B: 


B* b; 

D* d; 

b = new D ( ) ; 

d = (D*) b? // C-style cast 

d = dynami c_ca s t < D * > b; // C++ dynamic cast 

assert (d 1= NULL) ; 

The C++ standard only requires that dynamic casting be supported for classes 
with virtual functions. Dynamic casting is not supported by some older C++ 
compilers. See also R’lTl. 

eager evaluation Eager evaluation is the opposite of lazy evaluation. 

edge-centered field An edge-centered field is a field where there is a node associated 
with each edge of the mesh. 

Enterprise “Enterprise” is the name of the file format used for working with paged 
meshes and fields. See Chapter 23. 

exceptions Exceptions are a new feature of the C++ standard designed to support the 
handling of exceptional conditions. Exceptions are not supported by some older 
C++ compilers. 

face A cell c in a collection C is the face of a cell d 6 C if c can be defined in terms of 
a subset of the vertices of d. For example, a tetrahedron has triangle, edge, and 
vertex faces. See also proper face. 

facet k facet is a face. Typically facets refer to 2-cells, such as triangles and quadri- 
laterals. 

face-centered field A face-centered field is a field where a node is associated with 
each face in the mesh. Typically the faces refer to 2-cells, such as quadrilaterals 
and triangles. 
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factory A factory encapsulates the procedures required to build various types of ob- 
jects. For example, in FEL the class FEL_plot3d_f ield encapsulates the 
steps required to build each of the over 50 standard derived fields defined by 
PLOT3D. 

FAST The Flow Analysis Software Toolkit is a CFD visualization post-processing ap- 
plication developed at NASA Ames Research Center. 

FEL The Field Encapsulation Library . 

field A field represents a function within a domain via a finite set of nodes. Every field 
has a mesh that contains the location and organization of the nodes. 

friend class In C++, a class A can be declared to be a friend of another class B , en- 
abling A to have the same access rights to the member functions and data of B 
instances as B itself. 

general position The definition of general position is context dependent, i.e., depen- 
dent upon the particular application. Typically it is easier to define general po- 
sition in terms of what it is not: a configuration is not in general position if an 
infinitesimal perturbation can change how the configuration is classified. For 
example, 4 points in R 3 are not in general position if they are coplanar. 

grid The term grid is used as a synonym for mesh by many. For some, grid implies 
structured. Thus mesh is preferred to grid if one wishes to refer to unstructured 
objects, or both structured and unstructured objects. 

has-a relationship A has-a relationship denotes containment. Has-a is interchange- 
able with is-part-of or uses-for-implementation. See also is-a relationship. 

homeomorphic Two spaces are homeomorphic if they are topologically equivalent. 

hybrid mesh A hybrid mesh is a mesh that is partially structured and partially unstruc- 
tured. See Figure A. 1 . 



Figure A. 1 : An example of a hybrid mesh. 


IBLANK IBLANK is a standard from PLOT3D [WBPE92] for associating extra in- 
formation with each vertex in a mesh, using one integer per vertex. IBLANK val- 
ues can be used to provide hints about overlapping zones in multi-zone meshes, 
or to indicate that the node value associated with a vertex is not valid. 

incident cells Cells c and d are incident if c is a proper face of d, or vice versa. For 
example, there are vertices, edges, and quadrilaterals that are incident with a 
hexahedron in a mesh. 


instantiation Instantiation is the creation e "^ e generation at compile-Amk-lime 
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See also affine combination. 
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PL0T3D PL0T3D was originally a CFD visualization post processing tool developed 
by Pieter Buning of NASA Ames Research Center. PLOT3D lives on primarily 
in some of the standards it defined such as those for file formats and IBLANKs. 

proper face A cell c is a proper face of a cell d if c is a face of d and c ^ d. 

pure virtual function A pure virtual function of a class A is a member function that is 
required to be defined by subclasses of A. A class that has a pure virtual function 
is an abstract class. 

rectilinear mesh A rectilinear mesh is a structured mesh where the cells are aligned 
with the coordinate axes, but the spacing between adjacent vertices can be ir- 
regular. Rectilinear meshes are also known as irregular meshes or non-uniform 
meshes. See Figure A.2. 

reference counting Reference counting is a technique for tracking how many other 
objects are using (have a reference to) a particular object. Reference counting is 
typically part of a memory management strategy where objects are automatically 
deallocated (garbage collected) when their count goes to 0. Reference counting 
is also known as “use counting”. 

regular mesh A regular mesh is a structured mesh where the cells are aligned with 
the coordinate axes and the spacing between adjacent vertices is equal in each 
dimension. See Figure A.2. In some of the literature a regular mesh is considered 
to be a uniform mesh. 

RTTI Run-Time Type Identification is a new capability defined in the C++ standard 
allowing the user to query an object about its type. RTTI is required to be sup- 
ported only for classes which have virtual functions. RTTI is not supported by 

.... ^ some older C++ compilers. Some new C++ features, such as dynamic casting, 
5 depend upon Ike availability of RTTT. 

signature An operation’s signature is defined by its name, parameters, and return 
value. 

simplex A k-simplex is the convex hull of k + 7 affinely independent points. A 0- 
simplex is a vertex, a 1 -simplex is an edge, a 2-simplex is a triangle and a 3- 
simplex is a tetrahedron. 

simplicial decomposition Simplicial decomposition is the subdivision of mesh cells 
into simplices. For example, hexahedra would be subdivided into tetrahedra, 
and quadrilaterals would be subdivided into triangles. 

specialization In C++ templated functions, a specialization is an implementation for 
a specific template type that overrides the generic implementation provided by 
the template. 

star The star of a cell c in a cell complex C is the subset of C consisting of all the 
cells of which c is a face. 
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steady A steady mesh or field is one that does not change over time. See also unsteady. 

STL The Standard Template Library is a library of C++ templated classes supporting 
container data structures, such as sets. 

structured mesh In a structured mesh the vertices and cells are organized in a regular 
pattern. In one of the most typical types of structured meshes, the vertices can 
be addressed as if they were in a multi-dimensional array. In 3-d, all the 3- 
dimensional cells are hexahedral. Structured meshes can be further classified as 
uniform, regular, rectilinear and curvilinear. See Figure A.2. 



Uniform Regular Rectilinear 



Curvilinear 


Figure A.2: FEL structured mesh types. 


templates Templates are C++’s support for parameterized types. 

uniform mesh A uniform mesh is a structured mesh where the cells are aligned 
with the coordinate axes and the spacing between adjacent vertices is uniform 
throughout. See Figure A.2. 

unsteady An unsteady mesh or field is one that changes over time. See also steady. 

unstructured mesh An unstructured mesh in 3-d consists of polyhedral cells, with no 
implicit connectivity. The cells are not necessarily all tetrahedra, though this is 
one of the most common unstructured types. See Figure A.3. 



Figure A. 3: An unstructured mesh. 


use counting See reference counting. 

vertex -centered field A vertex-centered field is a field where a node is associated with 
each vertex in the mesh. The standard file formats defined by PLOT3D are all 
for vertex-centered fields. 
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virtual function A virtual function is a member function where the particular imple- 
mentation is selected at run-time, based on the type of the object for which the 
operation is called. C++ supports polymorphism via virtual functions. 

working set A working set is a subset of a much larger set, where maintaining a subset 
improves performance in some respect. For instance, the caching done in the 
memory hierarchies used by most microprocessor-based systems can be thought 
of as a performance improvement strategy based on working sets. 

zone A zone refers to a particular submesh in a multi-mesh. 
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